Linear Mixed Model were used.
Type: intact vs. scrambled)Category: faces vs. houses)Hemisphere: left vs. right)Duration: 17 vs. 200)## load libraries
library(tidyverse)
library(magrittr)
library(lme4)
library(lmerTest)
library(emmeans)
# options for calculating estimated marginal means
# lmer.df = c("kenward-roger", "satterthwaite", "asymptotic")
emm_options(lmer.df = "satterthwaite",
lmerTest.limit = 1e6,
opt.digits = FALSE) # , pbkrtest.limit =, disable.pbkrtest = FALSE
# general setting
saveCSV <- FALSE
ratio.count <- 0.8
RT.min <- 200
# set up the theme for plot and rainclound plot
plot_theme <- {
general_theme +
theme(legend.position = "right")
}
erp_theme <- {
general_theme +
theme(legend.position = "bottom")
}
ylimit.p1 <- c(2.5, -0.5)
ylimit.n170 <- c(1.5, -5.5) # y axis for plotting N170 amplitude
resp.labels = c("sure faces", "not sure faces", "guessing", "not sure houses", "sure houses")
Previous studies showed that the N170 amplitudes evoked by faces preseneted for longer durations were larger:
Duration <- c("17", "200")
Amplitudes <- c(-2.5, -5)
Hypotheses <- rep("means", 2)
df.lit <- data_frame(Duration, Amplitudes, Hypotheses)
## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.
plot.lit <- {
ggplot(data = df.lit, aes(y = Amplitudes, x = Duration, fill = Duration)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(color = "black") +
facet_grid(. ~ Hypotheses) +
scale_fill_manual(values = c(gray(.7), gray(.3))) +
scale_y_reverse(limits =c(0, -5.5), breaks = seq(0, -5, -1), expand= c(0, 0)) + # set the limit for y axis
geom_hline(yintercept = -2.5, linetype = "dashed", color = "red", size = 1) +
labs(x = "", y = expression(paste("Amplitudes (", mu, "V)"))) + # set the names for main, x and y axises , title = "Means",
# theme_bw() +
general_theme +
theme( # strip.text = element_text(color = "white"),
legend.position = "none",
axis.line = element_line(size = 3, colour = "black"))
}
plot.lit
There are two possible explainations for the aggregated mean amplitudes:
Duration <- rep(c(rep("17", 5), rep("200", 5) ), 2)
Amplitudes <- c(c(rep(-2.5, 5), rep(-5, 5)), c(rep(-.83, 3), rep(-5, 7)))
TrialNum <- rep(as.character(rep(c(1:5), 2)), 2)
Hypotheses <- c(rep("Possibility 1", 10), rep("Possibility 2", 10))
df.st <- data_frame(Duration, Amplitudes, TrialNum, Hypotheses)
plot.st <- {
ggplot(data = df.st, aes(y = Amplitudes, x = Duration, fill = Duration, color = TrialNum)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge") +
scale_fill_manual(values = c(gray(.7), gray(.3))) +
scale_color_manual(values = rep("black", 5)) +
facet_grid(. ~ Hypotheses) +
geom_hline(yintercept = -2.5, linetype = "dashed", color = "red", size = 1) +
scale_y_reverse(limits =c(0, -5.5), breaks = seq(0, -5, -1), expand= c(0, 0)) + # set the limit for y axis
labs(x = "Duration (ms)") + # set the names for main, x and y axises , title = "At the single-trial level:", , y = expression(paste("Amplitudes (", mu, "V)"))
theme_bw() +
general_theme +
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.line = element_line(size = 3, colour = "black"))
}
plot.st
raw.beha.205 <- {
"E205_35_Behavior.csv" %>%
file.path(folder.data, .) %>%
read_csv()
}
## Parsed with column specification:
## cols(
## .default = col_double(),
## ExperimentName = col_character(),
## Clock.Information = col_character(),
## Handedness = col_character(),
## ip_Address = col_character(),
## ResearcherID = col_character(),
## SessionDate = col_character(),
## SessionTime = col_time(format = ""),
## SessionTimeUtc = col_character(),
## Sex = col_character(),
## `FH[Block]` = col_character(),
## `Procedure[Block]` = col_character(),
## `Running[Block]` = col_character(),
## scrambleVar = col_character(),
## `stimCategory[Block]` = col_character(),
## `stimName[Block]` = col_character(),
## CellLabel = col_character(),
## `FH[Trial]` = col_character(),
## `Procedure[Trial]` = col_character(),
## `Running[Trial]` = col_character(),
## `stimCategory[Trial]` = col_character()
## # ... with 1 more columns
## )
## See spec(...) for full column specifications.
clean.beha.205 <- {
raw.beha.205 %>%
filter(`Procedure[Trial]` == "expProc") %>% # remove rows for practice and countdown
select(ExperimentName, SubjCode = Subject, blockID, Block, Trial, Age, Sex, Handedness, Type = `stimCategory[Trial]`, Category = `FH[Trial]`, Duration = `stimDuration[Trial]`, ACC = `Resp.ACC[Trial]`, Resp.RT = `Resp.RT[Trial]`, Resp = `Resp.RESP[Trial]`, Stimuli = `stimName[Trial]`) %>% # select and rename the columns (remove unuseful columns)
mutate(
SubjCode = factor(paste("P", SubjCode, sep = "")), # save SubjCode as factor
Type = substr(Type, 1, 1),
Type = recode_factor(Type, "N" = "intact", "S" = "scrambled"), # rename the levels of Type
Category = recode(Category, "hosue" = "house", "house" = "house", "face" = "face", .default = "face"), # rename "hosue" as "house"
Duration = factor(Duration), # save Duration as factor
Category = as.factor(Category),
Resp = as.factor(Resp),
blockID = as.factor(blockID)
) %>%
group_by(SubjCode, Type, Category, Duration, ACC) %>% # divide the data based on these conditions
mutate(
RT = Resp.RT + as.numeric(levels(Duration))[Duration],
RT.Z = scale(RT), # calculate the Z value for reaction times within each group
RT.Within3Z = ifelse(RT.Z <= 3 & RT.Z >= -3, 1, ifelse(RT.Z < -3 | RT.Z > 3, 0, NaN)) # if the Z values are between -3 and 3, will be marked as 1
) %>%
ungroup() # ungroup the data
}
# behavior.tidyup
head(clean.beha.205, 10)
# check the number of trials for 205
sum.count.205 <- {
clean.beha.205 %>%
group_by(SubjCode, Type, Category, Duration) %>%
summarize(Count_sum = n())
}
valid.count.205 <- {
clean.beha.205 %>%
filter(RT > RT.min) %>%
group_by(SubjCode, Type, Category, Duration) %>%
summarize(Count = n()) %>%
right_join(sum.count.205, by = c("SubjCode", "Type", "Category", "Duration")) %>%
mutate(ratio = Count / Count_sum) %>%
filter(ratio > ratio.count) %$%
SubjCode %>%
unique()
}
# remove the participants
clean.beha.205 %<>% filter(SubjCode %in% valid.count.205)
# remove participants whose performance in Key 1 was lower than 95%
valid.17 <- {
clean.beha.205 %>%
filter(Duration == 17, Resp == 1) %>% # Type == "intact",
group_by(SubjCode, Duration, Resp) %>% # Type,
summarize(Accuracy = mean(ACC), Count = n()) %>%
filter(Accuracy > 0.95) %$%
SubjCode
}
# clean.beha.205 %>%
# filter(Duration == 17, Resp == 1) %>% # Type == "intact",
# group_by(SubjCode, Duration, Resp) %>% # Type,
# summarize(Accuracy = mean(ACC), Count = n()) %>%
# filter(Accuracy < 0.95) %>%
# ungroup() %>%
# summarize(mean(Accuracy), min(Accuracy), max(Accuracy))
clean.beha.205 %<>% filter(SubjCode %in% valid.17)
subj.info.205 <- {
clean.beha.205 %>%
select(SubjCode, Age, Sex, Handedness) %>% # select the columns of participant inforamtion
distinct() # only keep the unique rows
}
subj.info.205
mean(subj.info.205$Age)
## [1] 22.32143
range(subj.info.205$Age)
## [1] 18 32
sum(subj.info.205$Sex == "female")
## [1] 20
tmp.count <- {
clean.beha.205 %>%
group_by(SubjCode, Type, Category, Duration) %>%
summarize(Count_sum = n())
}
avg.dist.long.205 <- {
clean.beha.205 %>%
# mutate(RespRateFace = if_else(Category == "face", ACC, if_else(Category == "house", as.numeric(!ACC), NaN))) %>%
group_by(SubjCode, Type, Category, Duration, Resp) %>%
summarize(Count = n()) %>%
right_join(tmp.count, by = c("SubjCode", "Type", "Category", "Duration")) %>%
mutate(RespRate = Count / Count_sum) %>%
ungroup()
}
# descriptive statistics of accuracy for plotting
desc.dist.205 <- {
avg.dist.long.205 %>%
group_by(Type, Category, Duration, Resp) %>%
summarize(Mean = mean(RespRate)
, N = n()
, SD = if_else(N > 1, sd(RespRate), 0)
, SE = SD/sqrt(N)
, SE.lo = Mean - SE, SE.hi = Mean + SE
, CI_lo = Mean - SE * 1.96
, CI_hi = Mean + SE * 1.96
, Median = median(RespRate)
, Lower = Mean - SD, Upper = Mean + SD # for rainclound plot
) %>%
ungroup() %>%
add_row(Type = "intact", Category = "house", Duration = 200, Resp = 2, Mean = 0, N = 0, SD = 0, SE = 0, SE.lo = 0,
SE.hi = 0, CI_lo = 0, CI_hi = 0, Median = 0, Lower = 0, Upper = 0)
}
avg.dist.wide.205 <- {
avg.dist.long.205 %>%
mutate(Conditions = paste(Type, Category, Duration, Resp, sep = "_")) %>%
select(SubjCode, Conditions, RespRate) %>%
spread(Conditions, RespRate)
}
knitr::kable(desc.dist.205, digits = 4)
| Type | Category | Duration | Resp | Mean | N | SD | SE | SE.lo | SE.hi | CI_lo | CI_hi | Median | Lower | Upper |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| intact | face | 17 | 1 | 0.4688 | 28 | 0.3306 | 0.0625 | 0.4063 | 0.5312 | 0.3463 | 0.5912 | 0.4469 | 0.1381 | 0.7994 |
| intact | face | 17 | 2 | 0.3205 | 28 | 0.2276 | 0.0430 | 0.2775 | 0.3636 | 0.2362 | 0.4049 | 0.3469 | 0.0929 | 0.5482 |
| intact | face | 17 | 3 | 0.1593 | 28 | 0.1631 | 0.0308 | 0.1284 | 0.1901 | 0.0988 | 0.2197 | 0.0969 | -0.0039 | 0.3224 |
| intact | face | 17 | 4 | 0.0478 | 28 | 0.0565 | 0.0107 | 0.0371 | 0.0584 | 0.0268 | 0.0687 | 0.0219 | -0.0087 | 0.1043 |
| intact | face | 17 | 5 | 0.0086 | 12 | 0.0129 | 0.0037 | 0.0049 | 0.0123 | 0.0013 | 0.0159 | 0.0031 | -0.0043 | 0.0215 |
| intact | face | 200 | 1 | 0.9768 | 28 | 0.0320 | 0.0061 | 0.9707 | 0.9828 | 0.9649 | 0.9887 | 0.9875 | 0.9447 | 1.0088 |
| intact | face | 200 | 2 | 0.0400 | 10 | 0.0403 | 0.0127 | 0.0273 | 0.0527 | 0.0150 | 0.0650 | 0.0188 | -0.0003 | 0.0803 |
| intact | face | 200 | 3 | 0.0208 | 9 | 0.0125 | 0.0042 | 0.0167 | 0.0250 | 0.0127 | 0.0290 | 0.0125 | 0.0083 | 0.0333 |
| intact | face | 200 | 4 | 0.0125 | 1 | 0.0000 | 0.0000 | 0.0125 | 0.0125 | 0.0125 | 0.0125 | 0.0125 | 0.0125 | 0.0125 |
| intact | face | 200 | 5 | 0.0125 | 4 | 0.0000 | 0.0000 | 0.0125 | 0.0125 | 0.0125 | 0.0125 | 0.0125 | 0.0125 | 0.0125 |
| intact | house | 17 | 1 | 0.0094 | 18 | 0.0086 | 0.0020 | 0.0073 | 0.0114 | 0.0054 | 0.0134 | 0.0062 | 0.0007 | 0.0180 |
| intact | house | 17 | 2 | 0.0550 | 28 | 0.0677 | 0.0128 | 0.0422 | 0.0678 | 0.0300 | 0.0801 | 0.0328 | -0.0126 | 0.1227 |
| intact | house | 17 | 3 | 0.3872 | 28 | 0.2582 | 0.0488 | 0.3384 | 0.4360 | 0.2915 | 0.4828 | 0.3234 | 0.1289 | 0.6454 |
| intact | house | 17 | 4 | 0.3988 | 28 | 0.2211 | 0.0418 | 0.3570 | 0.4406 | 0.3169 | 0.4807 | 0.4047 | 0.1776 | 0.6199 |
| intact | house | 17 | 5 | 0.1587 | 27 | 0.2304 | 0.0443 | 0.1143 | 0.2030 | 0.0718 | 0.2456 | 0.0438 | -0.0718 | 0.3891 |
| intact | house | 200 | 1 | 0.0188 | 6 | 0.0105 | 0.0043 | 0.0145 | 0.0230 | 0.0104 | 0.0271 | 0.0125 | 0.0083 | 0.0292 |
| intact | house | 200 | 3 | 0.0150 | 10 | 0.0079 | 0.0025 | 0.0125 | 0.0175 | 0.0101 | 0.0199 | 0.0125 | 0.0071 | 0.0229 |
| intact | house | 200 | 4 | 0.0331 | 20 | 0.0208 | 0.0046 | 0.0285 | 0.0378 | 0.0240 | 0.0422 | 0.0250 | 0.0123 | 0.0539 |
| intact | house | 200 | 5 | 0.9670 | 28 | 0.0262 | 0.0049 | 0.9620 | 0.9719 | 0.9573 | 0.9767 | 0.9750 | 0.9408 | 0.9931 |
| scrambled | face | 17 | 1 | 0.0156 | 12 | 0.0126 | 0.0036 | 0.0120 | 0.0193 | 0.0085 | 0.0228 | 0.0125 | 0.0030 | 0.0283 |
| scrambled | face | 17 | 2 | 0.0467 | 28 | 0.0538 | 0.0102 | 0.0365 | 0.0568 | 0.0267 | 0.0666 | 0.0375 | -0.0072 | 0.1005 |
| scrambled | face | 17 | 3 | 0.7143 | 28 | 0.2133 | 0.0403 | 0.6740 | 0.7546 | 0.6353 | 0.7933 | 0.7781 | 0.5010 | 0.9276 |
| scrambled | face | 17 | 4 | 0.2016 | 28 | 0.1637 | 0.0309 | 0.1706 | 0.2325 | 0.1409 | 0.2622 | 0.1656 | 0.0378 | 0.3653 |
| scrambled | face | 17 | 5 | 0.0616 | 14 | 0.1038 | 0.0277 | 0.0339 | 0.0893 | 0.0072 | 0.1160 | 0.0156 | -0.0422 | 0.1654 |
| scrambled | face | 200 | 1 | 0.0188 | 2 | 0.0088 | 0.0062 | 0.0125 | 0.0250 | 0.0065 | 0.0310 | 0.0188 | 0.0099 | 0.0276 |
| scrambled | face | 200 | 2 | 0.0292 | 12 | 0.0268 | 0.0077 | 0.0214 | 0.0369 | 0.0140 | 0.0443 | 0.0125 | 0.0023 | 0.0560 |
| scrambled | face | 200 | 3 | 0.9754 | 28 | 0.0332 | 0.0063 | 0.9692 | 0.9817 | 0.9631 | 0.9878 | 0.9875 | 0.9422 | 1.0087 |
| scrambled | face | 200 | 4 | 0.0216 | 11 | 0.0159 | 0.0048 | 0.0168 | 0.0264 | 0.0122 | 0.0310 | 0.0125 | 0.0057 | 0.0375 |
| scrambled | face | 200 | 5 | 0.0312 | 2 | 0.0265 | 0.0188 | 0.0125 | 0.0500 | -0.0055 | 0.0680 | 0.0312 | 0.0047 | 0.0578 |
| scrambled | house | 17 | 1 | 0.0144 | 10 | 0.0084 | 0.0026 | 0.0117 | 0.0170 | 0.0092 | 0.0196 | 0.0125 | 0.0060 | 0.0227 |
| scrambled | house | 17 | 2 | 0.0486 | 27 | 0.0520 | 0.0100 | 0.0386 | 0.0586 | 0.0290 | 0.0682 | 0.0375 | -0.0034 | 0.1006 |
| scrambled | house | 17 | 3 | 0.7002 | 28 | 0.2218 | 0.0419 | 0.6583 | 0.7421 | 0.6181 | 0.7824 | 0.7375 | 0.4785 | 0.9220 |
| scrambled | house | 17 | 4 | 0.2080 | 28 | 0.1729 | 0.0327 | 0.1754 | 0.2407 | 0.1440 | 0.2721 | 0.1531 | 0.0351 | 0.3809 |
| scrambled | house | 17 | 5 | 0.0695 | 16 | 0.1204 | 0.0301 | 0.0394 | 0.0996 | 0.0105 | 0.1285 | 0.0219 | -0.0509 | 0.1899 |
| scrambled | house | 200 | 1 | 0.0125 | 4 | 0.0000 | 0.0000 | 0.0125 | 0.0125 | 0.0125 | 0.0125 | 0.0125 | 0.0125 | 0.0125 |
| scrambled | house | 200 | 2 | 0.0250 | 11 | 0.0256 | 0.0077 | 0.0173 | 0.0327 | 0.0099 | 0.0401 | 0.0125 | -0.0006 | 0.0506 |
| scrambled | house | 200 | 3 | 0.9732 | 28 | 0.0345 | 0.0065 | 0.9667 | 0.9797 | 0.9604 | 0.9860 | 0.9750 | 0.9387 | 1.0077 |
| scrambled | house | 200 | 4 | 0.0240 | 12 | 0.0135 | 0.0039 | 0.0200 | 0.0279 | 0.0163 | 0.0316 | 0.0188 | 0.0104 | 0.0375 |
| scrambled | house | 200 | 5 | 0.0458 | 3 | 0.0473 | 0.0273 | 0.0185 | 0.0732 | -0.0077 | 0.0994 | 0.0250 | -0.0015 | 0.0932 |
| intact | house | 200 | 2 | 0.0000 | 0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
dist.RainPlot.205 <- {
ggplot(data = avg.dist.long.205, aes(y = RespRate, x = as.factor(Resp), fill = Category)) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .7) +
geom_point(aes(y = RespRate, color = Category), position = position_jitter(width = .15), size = .5, alpha = .8) +
geom_boxplot(aes(y = RespRate), width = .2, outlier.shape = NA, alpha = .7) +
geom_point(data = desc.dist.205, aes(y = Mean, x = Resp, color = Category), position = position_nudge(x = 0.3), size = 2.5) +
geom_errorbar(data = desc.dist.205, aes(y = Mean, ymin = Lower, ymax = Upper), position = position_nudge(x = 0.3), width = 0) +
facet_grid(Type ~ Duration) +
# facet_grid(. ~ Type, switch = "x") + # create two panels to show data
# scale_colour_grey() + # start = .1, end = .6, color for the contour
# scale_fill_grey() + # start = .3, end = .6, color for the fill
# scale_color_brewer(palette = "Set1") + # palette = "RdBu"
# scale_fill_brewer(palette = "Set1") + # palette = "RdBu"
labs(title = "Reaction times for 205", x = "Stimulus Type", y = "Reaction times (ms)", fill = "Duration(ms)", color = "Duration(ms)") + # set the names for main, x and y axises and the legend
# scale_x_discrete(labels=resp.labels) +
# coord_cartesian(ylim = c(0, 1.05)) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
theme_bw() + # the backgroud color
plot_theme
}
dist.RainPlot.205
dura.labels <- c("17ms", "200ms")
names(dura.labels) <- c(17, 200)
dist.ColuPlot.205 <- {
ggplot(data = desc.dist.205, aes(y = Mean, x = Resp, fill = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Type ~ Duration,
labeller = labeller(Duration = dura.labels)) +
geom_errorbar(mapping = aes(ymin = CI_lo, ymax = CI_hi), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
coord_cartesian(ylim = c(0,1.1)) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(title = "", x = "Responses", y = "Percentage", fill = "Category") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
erp_theme
}
dist.ColuPlot.205
df.ratio.205 <- {
clean.beha.205 %>%
mutate(isFace = if_else(Category == "face", 1, 0)) %>%
group_by(SubjCode, Type, Duration, Resp) %>%
summarize(ratio = mean(isFace), Count = n())
}
desc.ratio.205 <- {
df.ratio.205 %>%
group_by(Type, Duration, Resp) %>%
summarize(Mean = mean(ratio)
, N = n()
, SD = if_else(N > 1, sd(ratio), 0)
, SE = SD/sqrt(N)
, SE.lo = Mean - SE, SE.hi = Mean + SE
, CI_lo = Mean - SE * 1.96
, CI_hi = Mean + SE * 1.96
, Median = median(ratio)
, Lower = Mean - SD, Upper = Mean + SD # for rainclound plot
) %>%
ungroup()
}
df.ratio.205 %<>% sdif_coding_P203_resp
lmm.max.resp <- lmer(ratio ~ Type * Duration * Resp + (1 + Type_C + Dura_C + Type_Dura | SubjCode),
df.ratio.205)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -1.7e-02
print(summary(lmm.max.resp), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: ratio ~ Type * Duration * Resp + (1 + Type_C + Dura_C + Type_Dura | SubjCode)
## Data: df.ratio.205
##
## REML criterion at convergence: -88.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9487 -0.2868 -0.0109 0.1602 3.2327
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## SubjCode (Intercept) 0.0008608 0.02934
## Type_C 0.0009817 0.03133 1.00
## Dura_C 0.0046460 0.06816 0.72 0.72
## Type_Dura 0.0147012 0.12125 0.74 0.74 1.00
## Residual 0.0389778 0.19743
## Number of obs: 428, groups: SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.465725 0.012987 48.761546 35.860 < 2e-16 ***
## Type2-1 -0.009004 0.024204 94.303934 -0.372 0.710717
## Duration2-1 -0.020693 0.026829 43.357514 -0.771 0.444706
## Resp2-1 0.001016 0.038252 386.870545 0.027 0.978827
## Resp3-2 -0.281151 0.031701 366.105810 -8.869 < 2e-16 ***
## Resp4-3 -0.175603 0.029923 370.276934 -5.868 9.77e-09 ***
## Resp5-4 -0.062545 0.037613 381.720530 -1.663 0.097157 .
## Type2-1:Duration2-1 -0.154060 0.052316 49.742197 -2.945 0.004903 **
## Type2-1:Resp2-1 0.193941 0.076290 384.365785 2.542 0.011409 *
## Type2-1:Resp3-2 0.493453 0.063544 368.054343 7.765 8.15e-14 ***
## Type2-1:Resp4-3 0.272278 0.059862 370.772265 4.548 7.33e-06 ***
## Type2-1:Resp5-4 -0.087209 0.075176 380.506181 -1.160 0.246751
## Duration2-1:Resp2-1 0.285201 0.076504 386.781947 3.728 0.000222 ***
## Duration2-1:Resp3-2 -0.019693 0.063402 366.105810 -0.311 0.756283
## Duration2-1:Resp4-3 -0.197975 0.059846 370.276934 -3.308 0.001031 **
## Duration2-1:Resp5-4 -0.019521 0.075283 381.730772 -0.259 0.795542
## Type2-1:Duration2-1:Resp2-1 0.165273 0.152574 384.214140 1.083 0.279385
## Type2-1:Duration2-1:Resp3-2 -0.079215 0.127089 368.054342 -0.623 0.533471
## Type2-1:Duration2-1:Resp4-3 0.253131 0.119724 370.772265 2.114 0.035158 *
## Type2-1:Duration2-1:Resp5-4 -0.084537 0.150504 381.224642 -0.562 0.574656
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
lmm.zcp.resp <- update(lmm.max.resp,
formula = ratio ~ Type * Duration * Resp +
(1 + Type_C + Dura_C + Type_Dura || SubjCode))
## boundary (singular) fit: see ?isSingular
print(summary(lmm.zcp.resp), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: ratio ~ Type + Duration + Resp + (1 + Type_C + Dura_C + Type_Dura ||
## SubjCode) + Type:Duration + Type:Resp + Duration:Resp + Type:Duration:Resp
## Data: df.ratio.205
##
## REML criterion at convergence: -83.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9443 -0.2307 -0.0144 0.1260 3.3470
##
## Random effects:
## Groups Name Variance Std.Dev.
## SubjCode (Intercept) 1.719e-07 0.0004146
## SubjCode.1 Type_C 0.000e+00 0.0000000
## SubjCode.2 Dura_C 1.950e-03 0.0441571
## SubjCode.3 Type_Dura 4.430e-03 0.0665559
## Residual 4.057e-02 0.2014104
## Number of obs: 428, groups: SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.469614 0.011565 341.011588 40.607 < 2e-16 ***
## Type2-1 -0.001826 0.023125 394.027144 -0.079 0.937086
## Duration2-1 -0.013445 0.024611 44.049727 -0.546 0.587615
## Resp2-1 -0.004786 0.038271 386.579641 -0.125 0.900541
## Resp3-2 -0.280595 0.032260 374.465657 -8.698 < 2e-16 ***
## Resp4-3 -0.170501 0.030369 370.994699 -5.614 3.87e-08 ***
## Resp5-4 -0.057109 0.037598 380.493348 -1.519 0.129612
## Type2-1:Duration2-1 -0.140084 0.047982 44.817106 -2.919 0.005468 **
## Type2-1:Resp2-1 0.181034 0.076506 386.688023 2.366 0.018460 *
## Type2-1:Resp3-2 0.496577 0.064540 374.371374 7.694 1.28e-13 ***
## Type2-1:Resp4-3 0.280222 0.060739 370.985343 4.614 5.46e-06 ***
## Type2-1:Resp5-4 -0.075634 0.075188 380.497724 -1.006 0.315087
## Duration2-1:Resp2-1 0.276266 0.076546 386.822446 3.609 0.000348 ***
## Duration2-1:Resp3-2 -0.018579 0.064520 374.465655 -0.288 0.773545
## Duration2-1:Resp4-3 -0.187772 0.060738 370.994701 -3.092 0.002142 **
## Duration2-1:Resp5-4 -0.008628 0.075240 384.958509 -0.115 0.908761
## Type2-1:Duration2-1:Resp2-1 0.144795 0.153018 386.881798 0.946 0.344604
## Type2-1:Duration2-1:Resp3-2 -0.072967 0.129080 374.371376 -0.565 0.572219
## Type2-1:Duration2-1:Resp4-3 0.269019 0.121478 370.985343 2.215 0.027398 *
## Type2-1:Duration2-1:Resp5-4 -0.057947 0.150462 384.942732 -0.385 0.700356
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.zcp.resp))
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.3304 0.2192 0.002058 0
## Proportion of Variance 0.6943 0.3056 0.000030 0
## Cumulative Proportion 0.6943 1.0000 1.000000 1
lmm.rdc.resp <- update(lmm.zcp.resp,
formula = ratio ~ Type * Duration * Resp +
(0 + Dura_C + Type_Dura || SubjCode))
anova(lmm.rdc.resp, lmm.zcp.resp, refit = FALSE)
lmm.etd.resp <- update(lmm.zcp.resp,
formula = ratio ~ Type * Duration * Resp +
(0 + Dura_C + Type_Dura | SubjCode))
## boundary (singular) fit: see ?isSingular
summary(lmm.etd.resp)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: ratio ~ Type + Duration + Resp + (0 + Dura_C + Type_Dura | SubjCode) +
## Type:Duration + Type:Resp + Duration:Resp + Type:Duration:Resp
## Data: df.ratio.205
##
## REML criterion at convergence: -86
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8966 -0.2787 -0.0083 0.1667 3.3573
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## SubjCode Dura_C 0.003545 0.05954
## Type_Dura 0.011075 0.10524 1.00
## Residual 0.039931 0.19983
## Number of obs: 428, groups: SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.468898 0.011560 405.727681 40.564 < 2e-16 ***
## Type2-1 -0.002287 0.023116 405.312514 -0.099 0.921238
## Duration2-1 -0.014526 0.025787 48.991476 -0.563 0.575802
## Resp2-1 -0.004878 0.038157 401.292908 -0.128 0.898334
## Resp3-2 -0.280604 0.031988 386.549431 -8.772 < 2e-16 ***
## Resp4-3 -0.171915 0.030139 388.235257 -5.704 2.32e-08 ***
## Resp5-4 -0.056739 0.037493 399.754139 -1.513 0.130993
## Type2-1:Duration2-1 -0.140723 0.050486 56.564470 -2.787 0.007222 **
## Type2-1:Resp2-1 0.183236 0.076285 400.163399 2.402 0.016761 *
## Type2-1:Resp3-2 0.495156 0.063992 387.233395 7.738 8.88e-14 ***
## Type2-1:Resp4-3 0.277368 0.060278 388.268935 4.601 5.69e-06 ***
## Type2-1:Resp5-4 -0.075327 0.074980 399.423732 -1.005 0.315681
## Duration2-1:Resp2-1 0.272288 0.076320 401.315406 3.568 0.000404 ***
## Duration2-1:Resp3-2 -0.018597 0.063976 386.549431 -0.291 0.771451
## Duration2-1:Resp4-3 -0.190601 0.060278 388.235257 -3.162 0.001690 **
## Duration2-1:Resp5-4 -0.009927 0.075117 402.523833 -0.132 0.894923
## Type2-1:Duration2-1:Resp2-1 0.141611 0.152581 400.136396 0.928 0.353912
## Type2-1:Duration2-1:Resp3-2 -0.075808 0.127985 387.233395 -0.592 0.553982
## Type2-1:Duration2-1:Resp4-3 0.263312 0.120557 388.268935 2.184 0.029550 *
## Type2-1:Duration2-1:Resp5-4 -0.063515 0.150221 402.229698 -0.423 0.672661
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## boundary (singular) fit: see ?isSingular
lmm.etd1.resp <- update(lmm.etd.resp,
formula = ratio ~ Type * Duration * Resp +
(0 + Type_Dura | SubjCode))
summary(lmm.etd1.resp)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: ratio ~ Type + Duration + Resp + (0 + Type_Dura | SubjCode) +
## Type:Duration + Type:Resp + Duration:Resp + Type:Duration:Resp
## Data: df.ratio.205
##
## REML criterion at convergence: -83.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9667 -0.2385 -0.0109 0.1084 3.4064
##
## Random effects:
## Groups Name Variance Std.Dev.
## SubjCode Type_Dura 0.00390 0.06245
## Residual 0.04109 0.20270
## Number of obs: 428, groups: SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.469474 0.011599 400.842776 40.475 < 2e-16 ***
## Type2-1 -0.001081 0.023203 402.828429 -0.047 0.962865
## Duration2-1 -0.012959 0.023208 404.380853 -0.558 0.576907
## Resp2-1 -0.005719 0.038398 394.875793 -0.149 0.881683
## Resp3-2 -0.280367 0.032431 395.114720 -8.645 < 2e-16 ***
## Resp4-3 -0.169758 0.030528 390.645016 -5.561 4.99e-08 ***
## Resp5-4 -0.057450 0.037761 394.353731 -1.521 0.128954
## Type2-1:Duration2-1 -0.137427 0.047896 44.867934 -2.869 0.006253 **
## Type2-1:Resp2-1 0.181385 0.076846 401.877015 2.360 0.018734 *
## Type2-1:Resp3-2 0.496671 0.064835 389.636392 7.660 1.48e-13 ***
## Type2-1:Resp4-3 0.280445 0.061054 390.264928 4.593 5.89e-06 ***
## Type2-1:Resp5-4 -0.076454 0.075533 396.227120 -1.012 0.312064
## Duration2-1:Resp2-1 0.273275 0.076796 395.007070 3.558 0.000418 ***
## Duration2-1:Resp3-2 -0.018124 0.064862 395.114720 -0.279 0.780062
## Duration2-1:Resp4-3 -0.186286 0.061056 390.645016 -3.051 0.002436 **
## Duration2-1:Resp5-4 -0.006607 0.075536 396.632378 -0.087 0.930339
## Type2-1:Duration2-1:Resp2-1 0.143246 0.153695 402.003852 0.932 0.351887
## Type2-1:Duration2-1:Resp3-2 -0.072778 0.129671 389.636392 -0.561 0.574949
## Type2-1:Duration2-1:Resp4-3 0.269466 0.122109 390.264928 2.207 0.027912 *
## Type2-1:Duration2-1:Resp5-4 -0.056009 0.151094 398.314079 -0.371 0.711067
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
anova(lmm.etd1.resp, lmm.rdc.resp, refit = FALSE)
lmm.opt.resp <- lmm.etd1.resp
summary(lmm.opt.resp)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: ratio ~ Type + Duration + Resp + (0 + Type_Dura | SubjCode) +
## Type:Duration + Type:Resp + Duration:Resp + Type:Duration:Resp
## Data: df.ratio.205
##
## REML criterion at convergence: -83.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9667 -0.2385 -0.0109 0.1084 3.4064
##
## Random effects:
## Groups Name Variance Std.Dev.
## SubjCode Type_Dura 0.00390 0.06245
## Residual 0.04109 0.20270
## Number of obs: 428, groups: SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.469474 0.011599 400.842776 40.475 < 2e-16 ***
## Type2-1 -0.001081 0.023203 402.828429 -0.047 0.962865
## Duration2-1 -0.012959 0.023208 404.380853 -0.558 0.576907
## Resp2-1 -0.005719 0.038398 394.875793 -0.149 0.881683
## Resp3-2 -0.280367 0.032431 395.114720 -8.645 < 2e-16 ***
## Resp4-3 -0.169758 0.030528 390.645016 -5.561 4.99e-08 ***
## Resp5-4 -0.057450 0.037761 394.353731 -1.521 0.128954
## Type2-1:Duration2-1 -0.137427 0.047896 44.867934 -2.869 0.006253 **
## Type2-1:Resp2-1 0.181385 0.076846 401.877015 2.360 0.018734 *
## Type2-1:Resp3-2 0.496671 0.064835 389.636392 7.660 1.48e-13 ***
## Type2-1:Resp4-3 0.280445 0.061054 390.264928 4.593 5.89e-06 ***
## Type2-1:Resp5-4 -0.076454 0.075533 396.227120 -1.012 0.312064
## Duration2-1:Resp2-1 0.273275 0.076796 395.007070 3.558 0.000418 ***
## Duration2-1:Resp3-2 -0.018124 0.064862 395.114720 -0.279 0.780062
## Duration2-1:Resp4-3 -0.186286 0.061056 390.645016 -3.051 0.002436 **
## Duration2-1:Resp5-4 -0.006607 0.075536 396.632378 -0.087 0.930339
## Type2-1:Duration2-1:Resp2-1 0.143246 0.153695 402.003852 0.932 0.351887
## Type2-1:Duration2-1:Resp3-2 -0.072778 0.129671 389.636392 -0.561 0.574949
## Type2-1:Duration2-1:Resp4-3 0.269466 0.122109 390.264928 2.207 0.027912 *
## Type2-1:Duration2-1:Resp5-4 -0.056009 0.151094 398.314079 -0.371 0.711067
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
emm.resp <- emmeans(lmm.opt.resp, ~ Resp + Duration + Type)
as.data.frame(emm.resp) %>%
select(Type, Duration, Resp , everything())
resp.ColuPlot.205 <- {
ggplot(data = as.data.frame(emm.resp), aes(y = emmean, x = Resp)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Type ~ Duration,
labeller = labeller(Duration = dura.labels)) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(0.5, 1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
# coord_cartesian(ylim = c(0,1.12) ) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0), limits = c(-0.1, 1.15), breaks = seq(0, 1, 0.25)) + # remove the space between columns and x axis
labs(title = "", x = "Responses", y = "Percentage") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
plot_theme
}
resp.ColuPlot.205
contra.resp <- contrast(emmeans(lmm.opt.resp, ~ Duration | Resp + Type), "pairwise")
contra.resp[c(1, 5)]
## contrast Resp Type estimate SE df t.ratio p.value
## 17 - 200 1 intact -0.00760280 0.05449349 373.30 -0.140 0.8891
## 17 - 200 5 intact 0.07210303 0.05499259 374.51 1.311 0.1906
##
## Degrees-of-freedom method: satterthwaite
confint(contra.resp[c(1, 5)])
test(emm.resp, null = 0.5)
The raw mean amplitude of each trial for all participants were loaded.
# read the trial mean amplitude data
df.erp.E2 <- {
"E205_MeanAmp.csv" %>%
file.path(folder.data, .) %>%
read_csv() %>%
substr_Event() %>%
filter(SubjCode %in% valid.17) %>%
mutate(SubjCode = as.factor(SubjCode),
Hemisphere = as.factor(Hemisphere),
urResponse = as.factor(urResponse),
Stimuli = substr(stimName, 2, 6))
}
df.erp.E2 %>%
group_by(SubjCode, Stimuli, Hemisphere, Type, Category, Duration) %>%
summarize(count = n())
# dummy coding
df.erp.E2.all <- sdif_coding_P203_erp(df.erp.E2)
# only keep the data for amplitudes of the P1 (correct and incorrect erps)
df.erp.P1.E2 <- {
df.erp.E2.all %>%
filter(Component == "P1")
}
if (saveCSV) {
output.erp.P1.E2 <- file.path(folder.nesi, "E205_erp_P1.RData")
save(df.erp.P1.E2, file = output.erp.P1.E2)
}
# the structure of this dataset
head(df.erp.P1.E2, 10)
file.max.P1.E2 <- file.path(folder.nesi, "E205_P1_lmm_max.RData")
# fit the maximal model
if (!file.exists(file.max.P1.E2)) {
lmm.max.P1.E2 <- lmer(MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi +
Type_Cate_Dura_Hemi
| SubjCode) +
(1 + Type_C + Dura_C + Hemi_C +
Type_Dura + Type_Hemi + Dura_Hemi +
Type_Dura_Hemi
| Stimuli),
data = df.erp.P1.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e7)))
} else {
load(file.max.P1.E2)
}
print(summary(lmm.max.P1.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type * Category * Duration * Hemisphere + (1 + Type_C +
## Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura +
## Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi +
## Type_Dura_Hemi + Cate_Dura_Hemi + Type_Cate_Dura_Hemi | SubjCode) +
## (1 + Type_C + Dura_C + Hemi_C + Type_Dura + Type_Hemi + Dura_Hemi + Type_Dura_Hemi | Stimuli)
## Data: df.erp.P1.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 399168.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.4269 -0.5739 -0.0116 0.5710 9.2481
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.020242 0.14228
## Type_C 0.071284 0.26699 0.05
## Dura_C 0.072281 0.26885 0.40 -0.05
## Hemi_C 0.001384 0.03720 0.70 0.19 0.70
## Type_Dura 0.421342 0.64911 -0.10 0.67 0.12 -0.21
## Type_Hemi 0.011305 0.10632 -0.46 0.71 0.02 0.14 0.32
## Dura_Hemi 0.002092 0.04574 0.49 0.68 0.57 0.49 0.73 0.25
## Type_Dura_Hemi 0.075770 0.27526 -0.34 0.77 0.25 0.03 0.79 0.80 0.63
## SubjCode (Intercept) 1.797320 1.34064
## Type_C 0.113610 0.33706 0.10
## Cate_C 0.050124 0.22388 -0.36 -0.29
## Dura_C 0.321549 0.56705 0.44 -0.09 -0.36
## Hemi_C 0.594791 0.77123 0.31 0.11 0.48 -0.18
## Type_Cate 0.378197 0.61498 0.67 0.25 -0.36 0.38 0.10
## Type_Dura 0.336157 0.57979 0.05 0.22 0.05 0.29 0.02 0.41
## Cate_Dura 0.043316 0.20813 0.54 -0.40 0.09 0.18 0.55 0.41 -0.19
## Type_Hemi 0.034810 0.18657 0.02 0.47 -0.17 -0.02 -0.35 -0.27 -0.25 -0.45
## Cate_Hemi 0.011529 0.10737 0.00 0.55 0.23 -0.23 -0.06 0.03 0.17 -0.61 0.51
## Dura_Hemi 0.239585 0.48947 0.19 -0.08 -0.46 0.46 -0.37 0.23 0.19 0.27 0.22 -0.50
## Type_Cate_Dura 0.382269 0.61828 0.70 0.29 -0.67 0.53 -0.34 0.77 0.23 0.03 0.19 0.22 0.32
## Type_Cate_Hemi 0.127146 0.35658 0.00 -0.38 0.04 0.27 0.07 0.52 0.45 0.46 -0.83 -0.52 0.20 0.06
## Type_Dura_Hemi 0.133551 0.36545 -0.48 0.43 0.55 -0.36 0.34 -0.39 0.23 -0.26 0.30 0.29 -0.07 -0.59 -0.26
## Cate_Dura_Hemi 0.074409 0.27278 0.06 -0.44 0.18 0.21 -0.41 -0.16 -0.04 0.05 0.46 0.04 0.50 0.10 -0.18
## Type_Cate_Dura_Hemi 0.470022 0.68558 0.17 -0.64 0.44 0.34 -0.05 0.17 0.28 0.38 -0.16 -0.10 0.30 0.07 0.38
## Residual 15.678168 3.95957
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
## -0.03
## -0.12 0.75
##
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.262713 0.245906 29.244440 5.135 1.7e-05 ***
## Type2-1 -0.106661 0.076825 37.967608 -1.388 0.17312
## Category2-1 0.085615 0.062504 46.064875 1.370 0.17741
## Duration2-1 0.353467 0.113340 33.150333 3.119 0.00374 **
## Hemisphere2-1 0.117784 0.145144 29.038601 0.811 0.42368
## Type2-1:Category2-1 0.122340 0.145142 38.767656 0.843 0.40445
## Type2-1:Duration2-1 0.313789 0.146180 44.709815 2.147 0.03729 *
## Category2-1:Duration2-1 0.122743 0.099766 60.821395 1.230 0.22332
## Type2-1:Hemisphere2-1 0.011121 0.078696 72.888288 0.141 0.88801
## Category2-1:Hemisphere2-1 -0.041851 0.073108 178.077600 -0.572 0.56774
## Duration2-1:Hemisphere2-1 0.005918 0.113598 30.243456 0.052 0.95879
## Type2-1:Category2-1:Duration2-1 0.008769 0.231065 63.040636 0.038 0.96985
## Type2-1:Category2-1:Hemisphere2-1 0.190520 0.156105 79.923915 1.220 0.22588
## Type2-1:Duration2-1:Hemisphere2-1 -0.044562 0.158002 76.945281 -0.282 0.77867
## Category2-1:Duration2-1:Hemisphere2-1 -0.151587 0.148829 122.513631 -1.019 0.31043
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 0.123037 0.312596 84.251373 0.394 0.69487
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
file.zcp.P1.E2 <- file.path(folder.nesi, "E205_P1_lmm_zcp.RData")
# fit the zcp model
if (!file.exists(file.zcp.P1.E2)) {
lmm.zcp.P1.E2 <- lmer(MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi +
Type_Cate_Dura_Hemi
|| SubjCode) +
(1 + Type_C + Dura_C + Hemi_C +
Type_Dura + Type_Hemi + Dura_Hemi +
Type_Dura_Hemi
|| Stimuli),
data = df.erp.P1.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e7)))
} else {
load(file.zcp.P1.E2)
}
print(summary(lmm.zcp.P1.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type * Category * Duration * Hemisphere + (1 + Type_C +
## Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura +
## Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi +
## Type_Dura_Hemi + Cate_Dura_Hemi + Type_Cate_Dura_Hemi ||
## SubjCode) + (1 + Type_C + Dura_C + Hemi_C + Type_Dura + Type_Hemi + Dura_Hemi + Type_Dura_Hemi || Stimuli)
## Data: df.erp.P1.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 399290.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5143 -0.5726 -0.0123 0.5716 9.3014
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura_Hemi 1.932e-09 4.395e-05
## Stimuli.1 Dura_Hemi 2.685e-13 5.181e-07
## Stimuli.2 Type_Hemi 4.559e-13 6.752e-07
## Stimuli.3 Type_Dura 2.496e-01 4.996e-01
## Stimuli.4 Hemi_C 0.000e+00 0.000e+00
## Stimuli.5 Dura_C 5.021e-02 2.241e-01
## Stimuli.6 Type_C 3.631e-02 1.906e-01
## Stimuli.7 (Intercept) 1.472e-02 1.213e-01
## SubjCode Type_Cate_Dura_Hemi 3.519e-12 1.876e-06
## SubjCode.1 Cate_Dura_Hemi 1.916e-13 4.377e-07
## SubjCode.2 Type_Dura_Hemi 0.000e+00 0.000e+00
## SubjCode.3 Type_Cate_Hemi 0.000e+00 0.000e+00
## SubjCode.4 Type_Cate_Dura 0.000e+00 0.000e+00
## SubjCode.5 Dura_Hemi 2.525e-01 5.025e-01
## SubjCode.6 Cate_Hemi 0.000e+00 0.000e+00
## SubjCode.7 Type_Hemi 0.000e+00 0.000e+00
## SubjCode.8 Cate_Dura 3.003e-02 1.733e-01
## SubjCode.9 Type_Dura 3.098e-01 5.566e-01
## SubjCode.10 Type_Cate 2.762e-01 5.255e-01
## SubjCode.11 Hemi_C 6.411e-01 8.007e-01
## SubjCode.12 Dura_C 3.127e-01 5.592e-01
## SubjCode.13 Cate_C 4.097e-02 2.024e-01
## SubjCode.14 Type_C 1.055e-01 3.248e-01
## SubjCode.15 (Intercept) 1.781e+00 1.334e+00
## Residual 1.570e+01 3.963e+00
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.263e+00 2.446e-01 2.918e+01 5.162 1.59e-05 ***
## Type2-1 -1.068e-01 7.208e-02 3.472e+01 -1.482 0.1474
## Category2-1 8.635e-02 5.768e-02 4.827e+01 1.497 0.1409
## Duration2-1 3.532e-01 1.108e-01 3.198e+01 3.187 0.0032 **
## Hemisphere2-1 1.172e-01 1.503e-01 2.930e+01 0.780 0.4419
## Type2-1:Category2-1 1.227e-01 1.262e-01 4.055e+01 0.972 0.3366
## Type2-1:Duration2-1 3.127e-01 1.355e-01 4.016e+01 2.309 0.0262 *
## Category2-1:Duration2-1 1.223e-01 9.173e-02 4.585e+01 1.333 0.1891
## Type2-1:Hemisphere2-1 1.154e-02 6.999e-02 7.068e+04 0.165 0.8691
## Category2-1:Hemisphere2-1 -4.049e-02 6.999e-02 7.068e+04 -0.579 0.5629
## Duration2-1:Hemisphere2-1 5.761e-03 1.154e-01 3.009e+01 0.050 0.9605
## Type2-1:Category2-1:Duration2-1 7.014e-03 1.791e-01 9.725e+01 0.039 0.9688
## Type2-1:Category2-1:Hemisphere2-1 1.898e-01 1.400e-01 7.068e+04 1.356 0.1752
## Type2-1:Duration2-1:Hemisphere2-1 -4.251e-02 1.400e-01 7.068e+04 -0.304 0.7613
## Category2-1:Duration2-1:Hemisphere2-1 -1.505e-01 1.400e-01 7.068e+04 -1.075 0.2822
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 1.214e-01 2.800e-01 7.068e+04 0.434 0.6646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.zcp.P1.E2))
## $Stimuli
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## Standard deviation 0.1261 0.05655 0.04809 0.03062 1.109e-05 1.704e-07 1.308e-07 0
## Proportion of Variance 0.7114 0.14312 0.10350 0.04196 0.000e+00 0.000e+00 0.000e+00 0
## Cumulative Proportion 0.7114 0.85454 0.95804 1.00000 1.000e+00 1.000e+00 1.000e+00 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## Standard deviation 0.3368 0.2021 0.14111 0.14047 0.13262 0.12681 0.08195 0.05108 0.04373 4.734e-07 1.105e-07 0 0
## Proportion of Variance 0.4749 0.1710 0.08339 0.08264 0.07366 0.06734 0.02813 0.01093 0.00801 0.000e+00 0.000e+00 0 0
## Cumulative Proportion 0.4749 0.6459 0.72930 0.81194 0.88559 0.95294 0.98106 0.99199 1.00000 1.000e+00 1.000e+00 1 1
## [,14] [,15] [,16]
## Standard deviation 0 0 0
## Proportion of Variance 0 0 0
## Cumulative Proportion 1 1 1
file.rdc.P1.E2 <- file.path(folder.nesi, "E205_P1_lmm_rdc.RData")
# fit the rdc model
if (!file.exists(file.rdc.P1.E2)) {
# lmm.rdc.P1.E2.step <- step(lmm.rdc.P1.E2, reduce.fixed = FALSE)
lmm.rdc.P1.E2 <- update(lmm.zcp.P1.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Dura_Hemi
|| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
|| Stimuli))
} else {
load(file.rdc.P1.E2)
}
print(summary(lmm.rdc.P1.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C +
## Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura +
## Dura_Hemi || SubjCode) + (1 + Type_C + Dura_C + Type_Dura ||
## Stimuli) + Type:Category + Type:Duration + Category:Duration +
## Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere +
## Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere +
## Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.P1.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 399290.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5143 -0.5726 -0.0123 0.5716 9.3014
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura 0.24958 0.4996
## Stimuli.1 Dura_C 0.05021 0.2241
## Stimuli.2 Type_C 0.03631 0.1906
## Stimuli.3 (Intercept) 0.01472 0.1213
## SubjCode Dura_Hemi 0.25250 0.5025
## SubjCode.1 Cate_Dura 0.03003 0.1733
## SubjCode.2 Type_Dura 0.30984 0.5566
## SubjCode.3 Type_Cate 0.27617 0.5255
## SubjCode.4 Hemi_C 0.64112 0.8007
## SubjCode.5 Dura_C 0.31265 0.5592
## SubjCode.6 Cate_C 0.04097 0.2024
## SubjCode.7 Type_C 0.10546 0.3248
## SubjCode.8 (Intercept) 1.78060 1.3344
## Residual 15.70241 3.9626
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.263e+00 2.446e-01 2.918e+01 5.162 1.59e-05 ***
## Type2-1 -1.068e-01 7.208e-02 3.472e+01 -1.482 0.1474
## Category2-1 8.635e-02 5.768e-02 4.827e+01 1.497 0.1409
## Duration2-1 3.532e-01 1.108e-01 3.198e+01 3.187 0.0032 **
## Hemisphere2-1 1.172e-01 1.503e-01 2.930e+01 0.780 0.4419
## Type2-1:Category2-1 1.227e-01 1.262e-01 4.055e+01 0.972 0.3366
## Type2-1:Duration2-1 3.127e-01 1.355e-01 4.016e+01 2.309 0.0262 *
## Category2-1:Duration2-1 1.223e-01 9.173e-02 4.585e+01 1.333 0.1891
## Type2-1:Hemisphere2-1 1.154e-02 6.999e-02 7.068e+04 0.165 0.8691
## Category2-1:Hemisphere2-1 -4.049e-02 6.999e-02 7.068e+04 -0.579 0.5629
## Duration2-1:Hemisphere2-1 5.761e-03 1.154e-01 3.009e+01 0.050 0.9605
## Type2-1:Category2-1:Duration2-1 7.014e-03 1.791e-01 9.725e+01 0.039 0.9688
## Type2-1:Category2-1:Hemisphere2-1 1.898e-01 1.400e-01 7.068e+04 1.356 0.1752
## Type2-1:Duration2-1:Hemisphere2-1 -4.251e-02 1.400e-01 7.068e+04 -0.304 0.7613
## Category2-1:Duration2-1:Hemisphere2-1 -1.505e-01 1.400e-01 7.068e+04 -1.075 0.2822
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 1.214e-01 2.800e-01 7.068e+04 0.434 0.6646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.P1.E2, lmm.zcp.P1.E2, refit = FALSE)
file.etd.P1.E2 <- file.path(folder.nesi, "E205_P1_lmm_etd.RData")
# fit the etd model
if (!file.exists(file.etd.P1.E2)) {
lmm.etd.P1.E2 <- update(lmm.rdc.P1.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Dura_Hemi
| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
| Stimuli))
lmm.etd.P1.E2.2 <- re_fit(lmm.etd.P1.E2)
} else {
load(file.etd.P1.E2)
}
print(summary(lmm.etd.P1.E2.2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C +
## Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura +
## Dura_Hemi | SubjCode) + (1 + Type_C + Dura_C + Type_Dura |
## Stimuli) + Type:Category + Type:Duration + Category:Duration +
## Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere +
## Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere +
## Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.P1.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 399219.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.4381 -0.5727 -0.0121 0.5713 9.2907
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.02013 0.1419
## Type_C 0.06908 0.2628 0.06
## Dura_C 0.07107 0.2666 0.40 -0.06
## Type_Dura 0.41428 0.6436 -0.10 0.67 0.12
## SubjCode (Intercept) 1.79762 1.3408
## Type_C 0.11245 0.3353 0.11
## Cate_C 0.04079 0.2020 -0.35 -0.31
## Dura_C 0.32097 0.5665 0.44 -0.09 -0.37
## Hemi_C 0.61141 0.7819 0.30 0.11 0.50 -0.19
## Type_Cate 0.25555 0.5055 0.62 0.23 -0.23 0.31 0.22
## Type_Dura 0.33290 0.5770 0.05 0.22 0.07 0.29 0.03 0.43
## Cate_Dura 0.04222 0.2055 0.45 -0.45 0.13 0.10 0.58 0.41 -0.23
## Dura_Hemi 0.23709 0.4869 0.22 -0.08 -0.60 0.48 -0.43 0.20 0.16 0.14
## Residual 15.69460 3.9616
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.263e+00 2.459e-01 2.924e+01 5.135 1.7e-05 ***
## Type2-1 -1.065e-01 7.640e-02 3.752e+01 -1.394 0.17163
## Category2-1 8.568e-02 5.993e-02 4.572e+01 1.430 0.15958
## Duration2-1 3.540e-01 1.132e-01 3.306e+01 3.127 0.00367 **
## Hemisphere2-1 1.175e-01 1.470e-01 2.908e+01 0.799 0.43066
## Type2-1:Category2-1 1.216e-01 1.299e-01 4.603e+01 0.936 0.35391
## Type2-1:Duration2-1 3.135e-01 1.455e-01 4.433e+01 2.154 0.03669 *
## Category2-1:Duration2-1 1.222e-01 9.930e-02 6.054e+01 1.231 0.22312
## Type2-1:Hemisphere2-1 1.162e-02 6.997e-02 7.073e+04 0.166 0.86805
## Category2-1:Hemisphere2-1 -4.073e-02 6.997e-02 7.073e+04 -0.582 0.56048
## Duration2-1:Hemisphere2-1 6.368e-03 1.131e-01 2.998e+01 0.056 0.95549
## Type2-1:Category2-1:Duration2-1 6.142e-03 2.008e-01 7.769e+01 0.031 0.97567
## Type2-1:Category2-1:Hemisphere2-1 1.886e-01 1.399e-01 7.073e+04 1.348 0.17771
## Type2-1:Duration2-1:Hemisphere2-1 -4.223e-02 1.399e-01 7.073e+04 -0.302 0.76282
## Category2-1:Duration2-1:Hemisphere2-1 -1.508e-01 1.399e-01 7.073e+04 -1.078 0.28117
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 1.195e-01 2.799e-01 7.073e+04 0.427 0.66937
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.etd.P1.E2.2, lmm.rdc.P1.E2, refit = FALSE)
The reduced model (lmm.rdc.P1.E2) explained the data better than the extended model (lmm.rdc.P1.E2.2).
lmm.opt.P1.E205 <- lmm.rdc.P1.E2
# print(summary(lmm.opt.P1.E2), corr = FALSE)
summary(lmm.opt.P1.E205)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C +
## Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura +
## Dura_Hemi || SubjCode) + (1 + Type_C + Dura_C + Type_Dura ||
## Stimuli) + Type:Category + Type:Duration + Category:Duration +
## Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere +
## Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere +
## Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.P1.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 399290.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5143 -0.5726 -0.0123 0.5716 9.3014
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura 0.24958 0.4996
## Stimuli.1 Dura_C 0.05021 0.2241
## Stimuli.2 Type_C 0.03631 0.1906
## Stimuli.3 (Intercept) 0.01472 0.1213
## SubjCode Dura_Hemi 0.25250 0.5025
## SubjCode.1 Cate_Dura 0.03003 0.1733
## SubjCode.2 Type_Dura 0.30984 0.5566
## SubjCode.3 Type_Cate 0.27617 0.5255
## SubjCode.4 Hemi_C 0.64112 0.8007
## SubjCode.5 Dura_C 0.31265 0.5592
## SubjCode.6 Cate_C 0.04097 0.2024
## SubjCode.7 Type_C 0.10546 0.3248
## SubjCode.8 (Intercept) 1.78060 1.3344
## Residual 15.70241 3.9626
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.263e+00 2.446e-01 2.918e+01 5.162 1.59e-05 ***
## Type2-1 -1.068e-01 7.208e-02 3.472e+01 -1.482 0.1474
## Category2-1 8.635e-02 5.768e-02 4.827e+01 1.497 0.1409
## Duration2-1 3.532e-01 1.108e-01 3.198e+01 3.187 0.0032 **
## Hemisphere2-1 1.172e-01 1.503e-01 2.930e+01 0.780 0.4419
## Type2-1:Category2-1 1.227e-01 1.262e-01 4.055e+01 0.972 0.3366
## Type2-1:Duration2-1 3.127e-01 1.355e-01 4.016e+01 2.309 0.0262 *
## Category2-1:Duration2-1 1.223e-01 9.173e-02 4.585e+01 1.333 0.1891
## Type2-1:Hemisphere2-1 1.154e-02 6.999e-02 7.068e+04 0.165 0.8691
## Category2-1:Hemisphere2-1 -4.049e-02 6.999e-02 7.068e+04 -0.579 0.5629
## Duration2-1:Hemisphere2-1 5.761e-03 1.154e-01 3.009e+01 0.050 0.9605
## Type2-1:Category2-1:Duration2-1 7.014e-03 1.791e-01 9.725e+01 0.039 0.9688
## Type2-1:Category2-1:Hemisphere2-1 1.898e-01 1.400e-01 7.068e+04 1.356 0.1752
## Type2-1:Duration2-1:Hemisphere2-1 -4.251e-02 1.400e-01 7.068e+04 -0.304 0.7613
## Category2-1:Duration2-1:Hemisphere2-1 -1.505e-01 1.400e-01 7.068e+04 -1.075 0.2822
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 1.214e-01 2.800e-01 7.068e+04 0.434 0.6646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
# confint(lmm.opt.P1.E205, method="Wald")
# qqplot
qqplot_lmer(lmm.opt.P1.E205)
emm.P1.E205 <- emmeans(lmm.opt.P1.E205, ~ Type + Duration + Hemisphere + Category )
emm.P1.E205 %>%
as.data.frame() %>%
select(Type, Category, Duration, Hemisphere, everything())
# emmip(emm.P1.E205, Category ~ Duration | Type + Hemisphere, CIs = TRUE)
hemiLabel <- c("left", "right")
names(hemiLabel) <- c("Left", "Right")
p1.all.LinePlot = {
ggplot(data = data.frame(emm.P1.E205), aes(y = emmean, x = Duration, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(Type ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit.p1) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration (ms)", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# p1.all.LinePlot.slide <- {
# p1.all.LinePlot +
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('P1_all.png', p1.all.LinePlot.slide, width = 7, height = 4)
p1.all.LinePlot
contra.all.P1.gen <- contrast(emmeans(lmm.opt.P1.E205, ~ Duration | Hemisphere + Category + Type),
"pairwise", adjust = "Bonferroni")
contra.all.P1.gen[1:2]
## contrast Hemisphere Category Type estimate SE df t.ratio p.value
## 17 - 200 Left face intact -0.07107583 0.1643558 122.97 -0.432 1.0000
## 17 - 200 Right face intact -0.20370672 0.1643558 122.97 -1.239 0.4351
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 2 tests
confint(contra.all.P1.gen[1:2])
contra.all.P1.spec <- contrast(emmeans(lmm.opt.P1.E205, ~ Category + Duration | Hemisphere + Type),
interaction = "pairwise", adjust = "Bonferroni")
contra.all.P1.spec[1:2]
## Category_pairwise Duration_pairwise Hemisphere Type estimate SE df t.ratio p.value
## face - house 17 - 200 Left intact 0.22438787 0.1561692 266.23 1.437 0.3039
## face - house 17 - 200 Right intact 0.01316323 0.1561692 266.23 0.084 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 2 tests
confint(contra.all.P1.spec[1:2])
contra.all.P1.gen[1:8]
## contrast Hemisphere Category Type estimate SE df t.ratio p.value
## 17 - 200 Left face intact -0.0710758 0.1643558 122.97 -0.432 1.0000
## 17 - 200 Right face intact -0.2037067 0.1643558 122.97 -1.239 1.0000
## 17 - 200 Left house intact -0.2954637 0.1646024 123.70 -1.795 0.6007
## 17 - 200 Right house intact -0.2168700 0.1646024 123.70 -1.318 1.0000
## 17 - 200 Left face scrambled -0.4319184 0.1700147 140.79 -2.540 0.0972
## 17 - 200 Right face scrambled -0.4613435 0.1700147 140.79 -2.714 0.0599
## 17 - 200 Left house scrambled -0.6026303 0.1699668 140.64 -3.546 0.0043
## 17 - 200 Right house scrambled -0.5422116 0.1699668 140.64 -3.190 0.0140
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 8 tests
confint(contra.all.P1.gen[1:8])
contra.all.P1.spec
## Hemisphere = Left, Type = intact:
## Category_pairwise Duration_pairwise estimate SE df t.ratio p.value
## face - house 17 - 200 0.22438787 0.1561692 266.23 1.437 0.1519
##
## Hemisphere = Right, Type = intact:
## Category_pairwise Duration_pairwise estimate SE df t.ratio p.value
## face - house 17 - 200 0.01316323 0.1561692 266.23 0.084 0.9329
##
## Hemisphere = Left, Type = scrambled:
## Category_pairwise Duration_pairwise estimate SE df t.ratio p.value
## face - house 17 - 200 0.17071189 0.1675616 352.74 1.019 0.3090
##
## Hemisphere = Right, Type = scrambled:
## Category_pairwise Duration_pairwise estimate SE df t.ratio p.value
## face - house 17 - 200 0.08086813 0.1675616 352.74 0.483 0.6297
##
## Degrees-of-freedom method: satterthwaite
# confint(contra.all.P1.spec[1:4])
# only keep the data for amplitudes of the N170 (all erps)
df.erp.N170.E2 <- {
df.erp.E2.all %>%
filter(Component == "N170")
}
if (saveCSV) {
output.erp.N170.E2 <- file.path(folder.nesi, "E205_erp_N170.RData")
save(df.erp.N170.E2, file = output.erp.N170.E2)
}
# the structure of this dataset
head(df.erp.N170.E2, 10)
file.max.N170.E2 <- file.path(folder.nesi, "E205_N170_lmm_max.RData")
# fit the maximal model
if (!file.exists(file.max.N170.E2)) {
lmm.max.N170.E2 <- lmer(MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi +
Type_Cate_Dura_Hemi
| SubjCode) +
(1 + Type_C + Dura_C + Hemi_C +
Type_Dura + Type_Hemi + Dura_Hemi +
Type_Dura_Hemi
| Stimuli),
data = df.erp.N170.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e7)))
} else {
load(file.max.N170.E2)
}
print(summary(lmm.max.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type * Category * Duration * Hemisphere + (1 + Type_C +
## Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura +
## Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi +
## Type_Dura_Hemi + Cate_Dura_Hemi + Type_Cate_Dura_Hemi | SubjCode) +
## (1 + Type_C + Dura_C + Hemi_C + Type_Dura + Type_Hemi + Dura_Hemi + Type_Dura_Hemi | Stimuli)
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401341.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9374 -0.5879 -0.0046 0.5793 10.5966
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.017355 0.1317
## Type_C 0.099516 0.3155 0.23
## Dura_C 0.031794 0.1783 0.62 -0.33
## Hemi_C 0.004382 0.0662 0.32 0.31 -0.14
## Type_Dura 0.375660 0.6129 0.30 0.42 0.13 -0.61
## Type_Hemi 0.002172 0.0466 0.53 0.32 0.11 -0.40 0.89
## Dura_Hemi 0.010231 0.1011 0.49 0.63 0.04 0.87 -0.25 -0.19
## Type_Dura_Hemi 0.040130 0.2003 0.17 -0.51 0.38 -0.76 0.55 0.63 -0.77
## SubjCode (Intercept) 3.125355 1.7679
## Type_C 0.509356 0.7137 -0.23
## Cate_C 0.352785 0.5940 -0.27 0.64
## Dura_C 0.906309 0.9520 0.00 0.20 0.02
## Hemi_C 1.175660 1.0843 -0.21 -0.35 -0.02 -0.34
## Type_Cate 1.443237 1.2013 0.62 -0.57 -0.79 -0.01 -0.06
## Type_Dura 0.838114 0.9155 -0.04 0.59 0.36 -0.03 -0.53 -0.21
## Cate_Dura 0.441519 0.6645 0.21 0.22 0.53 -0.20 -0.16 -0.35 0.62
## Type_Hemi 0.520714 0.7216 -0.03 0.52 0.27 0.24 -0.59 -0.27 0.43 0.21
## Cate_Hemi 0.471528 0.6867 -0.02 0.07 -0.06 -0.01 -0.42 0.02 0.24 0.15 0.83
## Dura_Hemi 1.453583 1.2056 -0.15 0.03 -0.08 0.49 -0.18 0.01 0.17 0.01 0.31 0.26
## Type_Cate_Dura 1.523837 1.2344 0.39 -0.57 -0.70 -0.14 0.22 0.75 -0.62 -0.70 -0.49 -0.30 -0.14
## Type_Cate_Hemi 1.545688 1.2433 0.04 -0.20 0.00 -0.11 0.42 0.08 -0.10 -0.03 -0.86 -0.95 -0.19 0.30
## Type_Dura_Hemi 0.518336 0.7200 0.11 0.71 0.37 -0.08 -0.48 -0.15 0.60 0.21 0.65 0.41 -0.30 -0.34 -0.45
## Cate_Dura_Hemi 0.442789 0.6654 0.11 -0.17 -0.24 -0.43 -0.08 0.21 0.22 0.24 0.48 0.79 -0.08 -0.13 -0.66
## Type_Cate_Dura_Hemi 0.644599 0.8029 0.10 -0.36 0.05 0.29 0.24 0.00 -0.46 -0.11 -0.63 -0.70 -0.02 0.23 0.65
## Residual 16.105886 4.0132
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
## 0.37
## -0.66 -0.78
##
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.457986 0.323589 29.121560 -4.506 9.92e-05 ***
## Type2-1 1.035785 0.139571 32.818293 7.421 1.64e-08 ***
## Category2-1 0.953981 0.117833 32.584341 8.096 2.65e-09 ***
## Duration2-1 0.560394 0.178509 29.688145 3.139 0.00381 **
## Hemisphere2-1 0.002074 0.201247 29.015948 0.010 0.99185
## Type2-1:Category2-1 -1.932099 0.241068 34.283041 -8.015 2.29e-09 ***
## Type2-1:Duration2-1 1.014670 0.194070 36.684930 5.228 7.11e-06 ***
## Category2-1:Duration2-1 1.064759 0.146072 32.945035 7.289 2.32e-08 ***
## Type2-1:Hemisphere2-1 0.388565 0.149700 29.423292 2.596 0.01458 *
## Category2-1:Hemisphere2-1 0.188786 0.144784 30.757300 1.304 0.20194
## Duration2-1:Hemisphere2-1 -0.234009 0.231533 29.172466 -1.011 0.32047
## Type2-1:Category2-1:Duration2-1 -2.037947 0.299494 42.540100 -6.805 2.60e-08 ***
## Type2-1:Category2-1:Hemisphere2-1 -0.586785 0.267832 30.765442 -2.191 0.03616 *
## Type2-1:Duration2-1:Hemisphere2-1 0.454275 0.194635 38.826366 2.334 0.02487 *
## Category2-1:Duration2-1:Hemisphere2-1 0.046010 0.188081 41.085375 0.245 0.80796
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -0.529665 0.322344 53.344936 -1.643 0.10623
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
file.zcp.N170.E2 <- file.path(folder.nesi, "E205_N170_lmm_zcp.RData")
# fit the zcp model
if (!file.exists(file.zcp.N170.E2)) {
lmm.zcp.N170.E2 <- lmer(MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi + Cate_Dura_Hemi +
Type_Cate_Dura_Hemi
|| SubjCode) +
(1 + Type_C + Dura_C + Hemi_C +
Type_Dura + Type_Hemi + Dura_Hemi +
Type_Dura_Hemi
|| Stimuli),
data = df.erp.N170.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e7)))
} else {
load(file.zcp.N170.E2)
}
print(summary(lmm.zcp.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type * Category * Duration * Hemisphere + (1 + Type_C +
## Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura +
## Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi +
## Type_Dura_Hemi + Cate_Dura_Hemi + Type_Cate_Dura_Hemi ||
## SubjCode) + (1 + Type_C + Dura_C + Hemi_C + Type_Dura + Type_Hemi + Dura_Hemi + Type_Dura_Hemi || Stimuli)
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401644.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9367 -0.5852 -0.0038 0.5801 10.5465
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura_Hemi 1.275e-13 3.571e-07
## Stimuli.1 Dura_Hemi 2.294e-14 1.515e-07
## Stimuli.2 Type_Hemi 0.000e+00 0.000e+00
## Stimuli.3 Type_Dura 2.851e-01 5.340e-01
## Stimuli.4 Hemi_C 0.000e+00 0.000e+00
## Stimuli.5 Dura_C 8.721e-03 9.339e-02
## Stimuli.6 Type_C 6.560e-02 2.561e-01
## Stimuli.7 (Intercept) 1.123e-02 1.060e-01
## SubjCode Type_Cate_Dura_Hemi 1.072e-13 3.274e-07
## SubjCode.1 Cate_Dura_Hemi 3.997e-03 6.322e-02
## SubjCode.2 Type_Dura_Hemi 1.168e-01 3.417e-01
## SubjCode.3 Type_Cate_Hemi 1.076e+00 1.037e+00
## SubjCode.4 Type_Cate_Dura 1.287e+00 1.134e+00
## SubjCode.5 Dura_Hemi 1.510e+00 1.229e+00
## SubjCode.6 Cate_Hemi 2.930e-01 5.413e-01
## SubjCode.7 Type_Hemi 3.862e-01 6.215e-01
## SubjCode.8 Cate_Dura 4.000e-01 6.325e-01
## SubjCode.9 Type_Dura 7.758e-01 8.808e-01
## SubjCode.10 Type_Cate 1.289e+00 1.135e+00
## SubjCode.11 Hemi_C 1.186e+00 1.089e+00
## SubjCode.12 Dura_C 9.083e-01 9.530e-01
## SubjCode.13 Cate_C 3.241e-01 5.693e-01
## SubjCode.14 Type_C 4.811e-01 6.936e-01
## SubjCode.15 (Intercept) 3.128e+00 1.769e+00
## Residual 1.614e+01 4.017e+00
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.458e+00 3.236e-01 2.908e+01 -4.506 9.93e-05 ***
## Type2-1 1.034e+00 1.346e-01 3.173e+01 7.682 9.88e-09 ***
## Category2-1 9.531e-01 1.124e-01 3.175e+01 8.483 1.15e-09 ***
## Duration2-1 5.603e-01 1.779e-01 2.920e+01 3.150 0.00376 **
## Hemisphere2-1 1.826e-03 2.020e-01 2.904e+01 0.009 0.99285
## Type2-1:Category2-1 -1.930e+00 2.264e-01 3.295e+01 -8.524 7.62e-10 ***
## Type2-1:Duration2-1 1.012e+00 1.856e-01 3.524e+01 5.450 4.02e-06 ***
## Category2-1:Duration2-1 1.063e+00 1.371e-01 2.985e+01 7.748 1.25e-08 ***
## Type2-1:Hemisphere2-1 3.885e-01 1.338e-01 3.142e+01 2.903 0.00671 **
## Category2-1:Hemisphere2-1 1.895e-01 1.217e-01 3.280e+01 1.558 0.12891
## Duration2-1:Hemisphere2-1 -2.345e-01 2.353e-01 2.914e+01 -0.996 0.32721
## Type2-1:Category2-1:Duration2-1 -2.036e+00 2.780e-01 4.017e+01 -7.322 6.54e-09 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.890e-01 2.366e-01 3.446e+01 -2.489 0.01780 *
## Type2-1:Duration2-1:Hemisphere2-1 4.534e-01 1.550e-01 3.135e+01 2.925 0.00636 **
## Category2-1:Duration2-1:Hemisphere2-1 4.961e-02 1.424e-01 3.306e+01 0.348 0.72976
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.319e-01 2.838e-01 7.052e+04 -1.874 0.06095 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.zcp.N170.E2))
## $Stimuli
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## Standard deviation 0.1329 0.06376 0.02639 0.02325 8.89e-08 3.77e-08 0 0
## Proportion of Variance 0.7692 0.17698 0.03031 0.02353 0.00e+00 0.00e+00 0 0
## Cumulative Proportion 0.7692 0.94617 0.97647 1.00000 1.00e+00 1.00e+00 1 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## Standard deviation 0.4403 0.3059 0.28258 0.28236 0.27113 0.25820 0.23725 0.21927 0.17267 0.15744 0.15471 0.14172 0.13475
## Proportion of Variance 0.2376 0.1147 0.09788 0.09773 0.09011 0.08172 0.06899 0.05893 0.03655 0.03038 0.02934 0.02462 0.02226
## Cumulative Proportion 0.2376 0.3523 0.45020 0.54793 0.63804 0.71975 0.78875 0.84768 0.88423 0.91461 0.94395 0.96857 0.99083
## [,14] [,15] [,16]
## Standard deviation 0.08507 0.01574 8.15e-08
## Proportion of Variance 0.00887 0.00030 0.00e+00
## Cumulative Proportion 0.99970 1.00000 1.00e+00
file.rdc.N170.E2 <- file.path(folder.nesi, "E205_N170_lmm_rdc.RData")
# fit the rdc model
if (!file.exists(file.rdc.N170.E2)) {
# lmm.rdc.N170.E2.step <- step(lmm.rdc.N170.E2, reduce.fixed = FALSE)
lmm.rdc.N170.E2 <- update(lmm.zcp.N170.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi
|| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
|| Stimuli))
} else {
load(file.rdc.N170.E2)
}
print(summary(lmm.rdc.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C +
## Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura +
## Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi +
## Type_Dura_Hemi || SubjCode) + (1 + Type_C + Dura_C + Type_Dura ||
## Stimuli) + Type:Category + Type:Duration + Category:Duration +
## Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere +
## Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere +
## Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401644.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9367 -0.5852 -0.0037 0.5800 10.5465
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura 0.285125 0.53397
## Stimuli.1 Dura_C 0.008721 0.09339
## Stimuli.2 Type_C 0.065603 0.25613
## Stimuli.3 (Intercept) 0.011234 0.10599
## SubjCode Type_Dura_Hemi 0.116720 0.34164
## SubjCode.1 Type_Cate_Hemi 1.075684 1.03715
## SubjCode.2 Type_Cate_Dura 1.286463 1.13422
## SubjCode.3 Dura_Hemi 1.510497 1.22902
## SubjCode.4 Cate_Hemi 0.292453 0.54079
## SubjCode.5 Type_Hemi 0.386211 0.62146
## SubjCode.6 Cate_Dura 0.399993 0.63245
## SubjCode.7 Type_Dura 0.775818 0.88081
## SubjCode.8 Type_Cate 1.288519 1.13513
## SubjCode.9 Hemi_C 1.186257 1.08915
## SubjCode.10 Dura_C 0.908297 0.95305
## SubjCode.11 Cate_C 0.324088 0.56929
## SubjCode.12 Type_C 0.481105 0.69362
## SubjCode.13 (Intercept) 3.127930 1.76860
## Residual 16.137008 4.01709
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.458e+00 3.236e-01 2.908e+01 -4.506 9.93e-05 ***
## Type2-1 1.034e+00 1.346e-01 3.173e+01 7.682 9.88e-09 ***
## Category2-1 9.531e-01 1.124e-01 3.175e+01 8.483 1.15e-09 ***
## Duration2-1 5.603e-01 1.779e-01 2.920e+01 3.150 0.00376 **
## Hemisphere2-1 1.825e-03 2.020e-01 2.904e+01 0.009 0.99285
## Type2-1:Category2-1 -1.930e+00 2.264e-01 3.294e+01 -8.524 7.62e-10 ***
## Type2-1:Duration2-1 1.012e+00 1.856e-01 3.524e+01 5.450 4.02e-06 ***
## Category2-1:Duration2-1 1.063e+00 1.371e-01 2.985e+01 7.748 1.25e-08 ***
## Type2-1:Hemisphere2-1 3.885e-01 1.338e-01 3.143e+01 2.903 0.00671 **
## Category2-1:Hemisphere2-1 1.895e-01 1.216e-01 3.413e+01 1.559 0.12829
## Duration2-1:Hemisphere2-1 -2.345e-01 2.353e-01 2.914e+01 -0.996 0.32722
## Type2-1:Category2-1:Duration2-1 -2.036e+00 2.780e-01 4.017e+01 -7.322 6.54e-09 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.890e-01 2.366e-01 3.446e+01 -2.489 0.01780 *
## Type2-1:Duration2-1:Hemisphere2-1 4.534e-01 1.550e-01 3.135e+01 2.925 0.00636 **
## Category2-1:Duration2-1:Hemisphere2-1 4.963e-02 1.419e-01 7.054e+04 0.350 0.72655
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.319e-01 2.838e-01 7.054e+04 -1.874 0.06093 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.N170.E2, lmm.zcp.N170.E2, refit = FALSE)
file.etd.N170.E2 <- file.path(folder.nesi, "E205_N170_lmm_etd.RData")
# fit the etd model
if (!file.exists(file.etd.N170.E2)) {
lmm.etd.N170.E2 <- update(lmm.rdc.N170.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Cate_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Cate_Dura + Type_Hemi + Cate_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi + Type_Dura_Hemi
| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
| Stimuli))
} else {
load(file.etd.N170.E2)
}
print(summary(lmm.etd.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C +
## Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura +
## Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi +
## Type_Dura_Hemi | SubjCode) + (1 + Type_C + Dura_C + Type_Dura |
## Stimuli) + Type:Category + Type:Duration + Category:Duration +
## Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere +
## Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere +
## Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401372.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9468 -0.5869 -0.0042 0.5788 10.5930
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.01700 0.1304
## Type_C 0.09751 0.3123 0.22
## Dura_C 0.03059 0.1749 0.63 -0.36
## Type_Dura 0.36712 0.6059 0.32 0.43 0.12
## SubjCode (Intercept) 3.12552 1.7679
## Type_C 0.50943 0.7137 -0.23
## Cate_C 0.35290 0.5941 -0.27 0.64
## Dura_C 0.90612 0.9519 0.00 0.20 0.02
## Hemi_C 1.17515 1.0840 -0.20 -0.35 -0.02 -0.34
## Type_Cate 1.44334 1.2014 0.62 -0.57 -0.79 -0.01 -0.06
## Type_Dura 0.83973 0.9164 -0.04 0.59 0.36 -0.03 -0.53 -0.21
## Cate_Dura 0.44172 0.6646 0.21 0.21 0.53 -0.20 -0.16 -0.35 0.61
## Type_Hemi 0.51510 0.7177 -0.03 0.53 0.27 0.24 -0.60 -0.28 0.44 0.20
## Cate_Hemi 0.30258 0.5501 -0.05 0.12 0.00 0.12 -0.49 -0.04 0.21 0.11 0.86
## Dura_Hemi 1.45408 1.2059 -0.15 0.03 -0.08 0.49 -0.18 0.01 0.17 0.01 0.31 0.35
## Type_Cate_Dura 1.52322 1.2342 0.39 -0.58 -0.70 -0.14 0.22 0.75 -0.62 -0.70 -0.50 -0.33 -0.14
## Type_Cate_Hemi 1.16352 1.0787 0.04 -0.18 -0.02 -0.22 0.43 0.11 -0.02 0.01 -0.85 -0.95 -0.22 0.29
## Type_Dura_Hemi 0.49472 0.7034 0.11 0.72 0.39 -0.08 -0.49 -0.16 0.61 0.20 0.64 0.37 -0.31 -0.35 -0.37
## Residual 16.11624 4.0145
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.458e+00 3.236e-01 2.900e+01 -4.506 9.99e-05 ***
## Type2-1 1.036e+00 1.395e-01 2.389e+01 7.425 1.19e-07 ***
## Category2-1 9.539e-01 1.178e-01 2.705e+01 8.099 1.05e-08 ***
## Duration2-1 5.603e-01 1.785e-01 2.390e+01 3.140 0.00445 **
## Hemisphere2-1 1.499e-03 2.011e-01 1.289e+01 0.007 0.99417
## Type2-1:Category2-1 -1.932e+00 2.409e-01 2.846e+01 -8.021 8.72e-09 ***
## Type2-1:Duration2-1 1.015e+00 1.939e-01 3.399e+01 5.232 8.60e-06 ***
## Category2-1:Duration2-1 1.065e+00 1.459e-01 3.290e+01 7.298 2.29e-08 ***
## Type2-1:Hemisphere2-1 3.889e-01 1.490e-01 2.525e+01 2.610 0.01501 *
## Category2-1:Hemisphere2-1 1.901e-01 1.229e-01 3.634e+01 1.546 0.13075
## Duration2-1:Hemisphere2-1 -2.345e-01 2.313e-01 1.980e+01 -1.014 0.32287
## Type2-1:Category2-1:Duration2-1 -2.038e+00 2.988e-01 3.953e+01 -6.821 3.53e-08 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.894e-01 2.427e-01 3.582e+01 -2.429 0.02031 *
## Type2-1:Duration2-1:Hemisphere2-1 4.561e-01 1.913e-01 3.071e+01 2.384 0.02351 *
## Category2-1:Duration2-1:Hemisphere2-1 4.947e-02 1.418e-01 7.067e+04 0.349 0.72720
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.344e-01 2.836e-01 7.067e+04 -1.884 0.05956 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.etd.N170.E2))
## $Stimuli
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.1559 0.07239 0.04451 0.01604
## Proportion of Variance 0.7647 0.16489 0.06235 0.00809
## Cumulative Proportion 0.7647 0.92956 0.99191 1.00000
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## Standard deviation 0.5573 0.4522 0.3638 0.28544 0.23566 0.21812 0.16008 0.13781 0.06854 0.05289 0.0001626 2.121e-06 6.643e-08
## Proportion of Variance 0.3513 0.2313 0.1497 0.09216 0.06282 0.05382 0.02899 0.02148 0.00531 0.00316 0.0000000 0.000e+00 0.000e+00
## Cumulative Proportion 0.3513 0.5826 0.7322 0.82441 0.88724 0.94105 0.97004 0.99152 0.99684 1.00000 1.0000000 1.000e+00 1.000e+00
## [,14]
## Standard deviation 0
## Proportion of Variance 0
## Cumulative Proportion 1
file.etd1.N170.E2 <- file.path(folder.nesi, "E205_N170_lmm_etd1.RData")
# fit the etd1 model
if (!file.exists(file.etd1.N170.E2)) {
lmm.etd1.N170.E2 <- update(lmm.etd.N170.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Type_Hemi + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi
| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
| Stimuli))
} else {
load(file.etd1.N170.E2)
}
print(summary(lmm.etd1.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C +
## Dura_C + Hemi_C + Type_Cate + Type_Dura + Type_Hemi + Dura_Hemi +
## Type_Cate_Dura + Type_Cate_Hemi | SubjCode) + (1 + Type_C +
## Dura_C + Type_Dura | Stimuli) + Type:Category + Type:Duration +
## Category:Duration + Type:Hemisphere + Category:Hemisphere +
## Duration:Hemisphere + Type:Category:Duration + Type:Category:Hemisphere +
## Type:Duration:Hemisphere + Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401774.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9501 -0.5858 -0.0039 0.5816 10.6037
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.01651 0.1285
## Type_C 0.09623 0.3102 0.21
## Dura_C 0.03051 0.1747 0.64 -0.36
## Type_Dura 0.36483 0.6040 0.33 0.42 0.13
## SubjCode (Intercept) 3.12440 1.7676
## Type_C 0.50975 0.7140 -0.23
## Dura_C 0.90542 0.9515 0.00 0.20
## Hemi_C 1.15652 1.0754 -0.21 -0.34 -0.34
## Type_Cate 1.75252 1.3238 0.62 -0.59 -0.03 -0.07
## Type_Dura 0.84363 0.9185 -0.04 0.60 -0.03 -0.53 -0.20
## Type_Hemi 0.38161 0.6177 -0.06 0.42 0.30 -0.56 -0.28 0.34
## Dura_Hemi 1.48018 1.2166 -0.15 0.01 0.49 -0.17 0.03 0.15 0.43
## Type_Cate_Dura 1.31094 1.1450 0.30 -0.44 -0.10 0.25 0.55 -0.66 -0.49 -0.17
## Type_Cate_Hemi 1.69344 1.3013 0.04 -0.18 -0.21 0.46 0.09 -0.05 -0.90 -0.25 0.33
## Residual 16.22266 4.0277
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.457e+00 3.235e-01 2.912e+01 -4.504 9.97e-05 ***
## Type2-1 1.035e+00 1.395e-01 3.271e+01 7.417 1.69e-08 ***
## Category2-1 9.523e-01 4.574e-02 7.773e+01 20.822 < 2e-16 ***
## Duration2-1 5.637e-01 1.784e-01 2.968e+01 3.160 0.00362 **
## Hemisphere2-1 2.036e-03 1.995e-01 2.905e+01 0.010 0.99193
## Type2-1:Category2-1 -1.928e+00 2.613e-01 3.344e+01 -7.379 1.64e-08 ***
## Type2-1:Duration2-1 1.017e+00 1.943e-01 3.637e+01 5.235 7.13e-06 ***
## Category2-1:Duration2-1 1.060e+00 8.118e-02 7.754e+01 13.057 < 2e-16 ***
## Type2-1:Hemisphere2-1 3.891e-01 1.333e-01 3.398e+01 2.918 0.00620 **
## Category2-1:Hemisphere2-1 1.909e-01 7.114e-02 7.070e+04 2.683 0.00730 **
## Duration2-1:Hemisphere2-1 -2.350e-01 2.332e-01 2.916e+01 -1.007 0.32199
## Type2-1:Category2-1:Duration2-1 -2.031e+00 2.867e-01 4.327e+01 -7.086 9.37e-09 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.876e-01 2.769e-01 3.445e+01 -2.122 0.04111 *
## Type2-1:Duration2-1:Hemisphere2-1 4.548e-01 1.423e-01 7.070e+04 3.196 0.00139 **
## Category2-1:Duration2-1:Hemisphere2-1 4.671e-02 1.423e-01 7.070e+04 0.328 0.74268
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.362e-01 2.846e-01 7.070e+04 -1.884 0.05956 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.etd1.N170.E2))
## $Stimuli
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.1548 0.07198 0.04396 0.01573
## Proportion of Variance 0.7650 0.16543 0.06171 0.00790
## Cumulative Proportion 0.7650 0.93039 0.99210 1.00000
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## Standard deviation 0.5368 0.4428 0.3314 0.28204 0.22279 0.21459 0.1596 0.12076 0.04202 6.325e-08
## Proportion of Variance 0.3553 0.2418 0.1354 0.09807 0.06119 0.05677 0.0314 0.01798 0.00218 0.000e+00
## Cumulative Proportion 0.3553 0.5970 0.7324 0.83048 0.89168 0.94845 0.9798 0.99782 1.00000 1.000e+00
file.etd2.N170.E2 <- file.path(folder.nesi, "E205_N170_lmm_etd2.RData")
# fit the etd2 model
if (!file.exists(file.etd2.N170.E2)) {
lmm.etd2.N170.E2 <- update(lmm.etd1.N170.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Type_C + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi
| SubjCode) +
(1 + Type_C + Dura_C +
Type_Dura
| Stimuli))
} else {
load(file.etd2.N170.E2)
}
print(summary(lmm.etd2.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C +
## Dura_C + Hemi_C + Type_Cate + Type_Dura + Dura_Hemi + Type_Cate_Dura +
## Type_Cate_Hemi | SubjCode) + (1 + Type_C + Dura_C + Type_Dura |
## Stimuli) + Type:Category + Type:Duration + Category:Duration +
## Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere +
## Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere +
## Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401871.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9303 -0.5860 -0.0048 0.5817 10.5986
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.002805 0.05297
## Type_C 0.092295 0.30380 1.00
## Dura_C 0.011615 0.10777 -0.32 -0.32
## Type_Dura 0.399586 0.63213 0.44 0.44 0.08
## SubjCode (Intercept) 3.123796 1.76743
## Type_C 0.509423 0.71374 -0.23
## Dura_C 0.905459 0.95156 0.00 0.20
## Hemi_C 1.218893 1.10403 -0.20 -0.35 -0.35
## Type_Cate 1.753881 1.32434 0.62 -0.59 -0.03 -0.05
## Type_Dura 0.838924 0.91593 -0.04 0.60 -0.03 -0.53 -0.20
## Dura_Hemi 1.601393 1.26546 -0.15 0.04 0.49 -0.23 0.00 0.17
## Type_Cate_Dura 1.313151 1.14593 0.30 -0.44 -0.10 0.27 0.55 -0.67 -0.21
## Type_Cate_Hemi 1.594705 1.26282 0.04 -0.18 -0.21 0.50 0.09 -0.04 -0.31 0.36
## Residual 16.256515 4.03194
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.457e+00 3.232e-01 2.902e+01 -4.507 9.93e-05 ***
## Type2-1 1.034e+00 1.393e-01 3.248e+01 7.425 1.73e-08 ***
## Category2-1 9.525e-01 3.754e-02 4.821e+01 25.376 < 2e-16 ***
## Duration2-1 5.641e-01 1.778e-01 2.926e+01 3.174 0.00353 **
## Hemisphere2-1 2.941e-03 2.047e-01 2.905e+01 0.014 0.98864
## Type2-1:Category2-1 -1.929e+00 2.611e-01 3.316e+01 -7.387 1.69e-08 ***
## Type2-1:Duration2-1 1.016e+00 1.950e-01 3.715e+01 5.208 7.32e-06 ***
## Category2-1:Duration2-1 1.060e+00 7.520e-02 1.147e+02 14.100 < 2e-16 ***
## Type2-1:Hemisphere2-1 3.885e-01 7.122e-02 7.080e+04 5.455 4.92e-08 ***
## Category2-1:Hemisphere2-1 1.916e-01 7.122e-02 7.080e+04 2.690 0.00714 **
## Duration2-1:Hemisphere2-1 -2.349e-01 2.418e-01 2.922e+01 -0.972 0.33917
## Type2-1:Category2-1:Duration2-1 -2.031e+00 2.899e-01 4.504e+01 -7.006 9.93e-09 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.857e-01 2.710e-01 3.421e+01 -2.161 0.03776 *
## Type2-1:Duration2-1:Hemisphere2-1 4.493e-01 1.424e-01 7.080e+04 3.155 0.00161 **
## Category2-1:Duration2-1:Hemisphere2-1 4.771e-02 1.424e-01 7.080e+04 0.335 0.73767
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.345e-01 2.849e-01 7.080e+04 -1.876 0.06065 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.etd2.N170.E2))
## $Stimuli
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.1611 0.0678 0.02407 1.847e-18
## Proportion of Variance 0.8338 0.1476 0.01861 0.000e+00
## Cumulative Proportion 0.8338 0.9814 1.00000 1.000e+00
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## Standard deviation 0.5344 0.4352 0.3311 0.26908 0.2225 0.21374 0.15795 0.11781 2.884e-05
## Proportion of Variance 0.3610 0.2394 0.1386 0.09153 0.0626 0.05775 0.03154 0.01755 0.000e+00
## Cumulative Proportion 0.3610 0.6004 0.7390 0.83057 0.8932 0.95092 0.98245 1.00000 1.000e+00
file.etd3.N170.E2 <- file.path(folder.nesi, "E205_N170_lmm_etd3.RData")
# fit the etd3 model
if (!file.exists(file.etd3.N170.E2)) {
lmm.etd3.N170.E2 <- update(lmm.etd2.N170.E2,
formula = MeanAmp ~ Type * Category * Duration * Hemisphere +
(1 + Dura_C + Hemi_C +
Type_Cate + Type_Dura + Dura_Hemi +
Type_Cate_Dura + Type_Cate_Hemi
| SubjCode) +
(0 + Type_C + Dura_C +
Type_Dura
| Stimuli))
} else {
load(file.etd3.N170.E2)
}
print(summary(lmm.etd3.N170.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Dura_C +
## Hemi_C + Type_Cate + Type_Dura + Dura_Hemi + Type_Cate_Dura +
## Type_Cate_Hemi | SubjCode) + (0 + Type_C + Dura_C + Type_Dura |
## Stimuli) + Type:Category + Type:Duration + Category:Duration +
## Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere +
## Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere +
## Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 402199
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9564 -0.5894 -0.0039 0.5825 10.5129
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli Type_C 0.09270 0.3045
## Dura_C 0.01805 0.1343 -0.66
## Type_Dura 0.38023 0.6166 0.41 -0.14
## SubjCode (Intercept) 3.14852 1.7744
## Dura_C 0.92884 0.9638 -0.02
## Hemi_C 1.21852 1.1039 -0.19 -0.37
## Type_Cate 1.75483 1.3247 0.63 -0.07 -0.05
## Type_Dura 0.54862 0.7407 0.15 -0.23 -0.34 0.29
## Dura_Hemi 1.60033 1.2650 -0.15 0.49 -0.23 0.00 0.17
## Type_Cate_Dura 1.31892 1.1484 0.30 -0.13 0.27 0.55 -0.43 -0.20
## Type_Cate_Hemi 1.59093 1.2613 0.05 -0.22 0.50 0.09 0.12 -0.31 0.36
## Residual 16.34600 4.0430
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.456e+00 3.245e-01 2.900e+01 -4.488 0.000105 ***
## Type2-1 1.031e+00 4.934e-02 7.879e+01 20.905 < 2e-16 ***
## Category2-1 9.524e-01 3.571e-02 7.090e+04 26.665 < 2e-16 ***
## Duration2-1 5.661e-01 1.802e-01 2.942e+01 3.142 0.003812 **
## Hemisphere2-1 2.973e-03 2.047e-01 2.905e+01 0.015 0.988510
## Type2-1:Category2-1 -1.924e+00 2.612e-01 3.328e+01 -7.367 1.76e-08 ***
## Type2-1:Duration2-1 1.011e+00 1.678e-01 4.324e+01 6.025 3.30e-07 ***
## Category2-1:Duration2-1 1.060e+00 7.749e-02 1.154e+02 13.685 < 2e-16 ***
## Type2-1:Hemisphere2-1 3.885e-01 7.142e-02 7.083e+04 5.439 5.36e-08 ***
## Category2-1:Hemisphere2-1 1.916e-01 7.141e-02 7.083e+04 2.683 0.007308 **
## Duration2-1:Hemisphere2-1 -2.349e-01 2.418e-01 2.921e+01 -0.972 0.339161
## Type2-1:Category2-1:Duration2-1 -2.021e+00 2.888e-01 4.417e+01 -7.000 1.12e-08 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.858e-01 2.710e-01 3.424e+01 -2.162 0.037722 *
## Type2-1:Duration2-1:Hemisphere2-1 4.492e-01 1.428e-01 7.083e+04 3.145 0.001661 **
## Category2-1:Duration2-1:Hemisphere2-1 4.772e-02 1.428e-01 7.083e+04 0.334 0.738281
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.343e-01 2.857e-01 7.083e+04 -1.870 0.061423 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.etd3.N170.E2, lmm.rdc.N170.E2, refit = FALSE)
lmm.opt.N170.E205 <- lmm.rdc.N170.E2
print(summary(lmm.opt.N170.E205), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Type + Category + Duration + Hemisphere + (1 + Type_C +
## Cate_C + Dura_C + Hemi_C + Type_Cate + Type_Dura + Cate_Dura +
## Type_Hemi + Cate_Hemi + Dura_Hemi + Type_Cate_Dura + Type_Cate_Hemi +
## Type_Dura_Hemi || SubjCode) + (1 + Type_C + Dura_C + Type_Dura ||
## Stimuli) + Type:Category + Type:Duration + Category:Duration +
## Type:Hemisphere + Category:Hemisphere + Duration:Hemisphere +
## Type:Category:Duration + Type:Category:Hemisphere + Type:Duration:Hemisphere +
## Category:Duration:Hemisphere + Type:Category:Duration:Hemisphere
## Data: df.erp.N170.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+07))
##
## REML criterion at convergence: 401644.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9367 -0.5852 -0.0037 0.5800 10.5465
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Type_Dura 0.285125 0.53397
## Stimuli.1 Dura_C 0.008721 0.09339
## Stimuli.2 Type_C 0.065603 0.25613
## Stimuli.3 (Intercept) 0.011234 0.10599
## SubjCode Type_Dura_Hemi 0.116720 0.34164
## SubjCode.1 Type_Cate_Hemi 1.075684 1.03715
## SubjCode.2 Type_Cate_Dura 1.286463 1.13422
## SubjCode.3 Dura_Hemi 1.510497 1.22902
## SubjCode.4 Cate_Hemi 0.292453 0.54079
## SubjCode.5 Type_Hemi 0.386211 0.62146
## SubjCode.6 Cate_Dura 0.399993 0.63245
## SubjCode.7 Type_Dura 0.775818 0.88081
## SubjCode.8 Type_Cate 1.288519 1.13513
## SubjCode.9 Hemi_C 1.186257 1.08915
## SubjCode.10 Dura_C 0.908297 0.95305
## SubjCode.11 Cate_C 0.324088 0.56929
## SubjCode.12 Type_C 0.481105 0.69362
## SubjCode.13 (Intercept) 3.127930 1.76860
## Residual 16.137008 4.01709
## Number of obs: 71282, groups: Stimuli, 80; SubjCode, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.458e+00 3.236e-01 2.908e+01 -4.506 9.93e-05 ***
## Type2-1 1.034e+00 1.346e-01 3.173e+01 7.682 9.88e-09 ***
## Category2-1 9.531e-01 1.124e-01 3.175e+01 8.483 1.15e-09 ***
## Duration2-1 5.603e-01 1.779e-01 2.920e+01 3.150 0.00376 **
## Hemisphere2-1 1.825e-03 2.020e-01 2.904e+01 0.009 0.99285
## Type2-1:Category2-1 -1.930e+00 2.264e-01 3.294e+01 -8.524 7.62e-10 ***
## Type2-1:Duration2-1 1.012e+00 1.856e-01 3.524e+01 5.450 4.02e-06 ***
## Category2-1:Duration2-1 1.063e+00 1.371e-01 2.985e+01 7.748 1.25e-08 ***
## Type2-1:Hemisphere2-1 3.885e-01 1.338e-01 3.143e+01 2.903 0.00671 **
## Category2-1:Hemisphere2-1 1.895e-01 1.216e-01 3.413e+01 1.559 0.12829
## Duration2-1:Hemisphere2-1 -2.345e-01 2.353e-01 2.914e+01 -0.996 0.32722
## Type2-1:Category2-1:Duration2-1 -2.036e+00 2.780e-01 4.017e+01 -7.322 6.54e-09 ***
## Type2-1:Category2-1:Hemisphere2-1 -5.890e-01 2.366e-01 3.446e+01 -2.489 0.01780 *
## Type2-1:Duration2-1:Hemisphere2-1 4.534e-01 1.550e-01 3.135e+01 2.925 0.00636 **
## Category2-1:Duration2-1:Hemisphere2-1 4.963e-02 1.419e-01 7.054e+04 0.350 0.72655
## Type2-1:Category2-1:Duration2-1:Hemisphere2-1 -5.319e-01 2.838e-01 7.054e+04 -1.874 0.06093 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# qqplot
qqplot_lmer(lmm.opt.N170.E205)
emm.N170.E205 <- emmeans(lmm.opt.N170.E205, ~ Duration + Category + Hemisphere + Type)
emm.N170.E205 %>%
as.data.frame() %>%
select(Type, Category, Duration, Hemisphere, everything())
# emmip(emm.N170.E205, Category ~ Duration | Type + Hemisphere, CIs = TRUE) +
# scale_y_continuous(trans = "reverse")
n170.all.LinePlot = {
ggplot(data = data.frame(emm.N170.E205), aes(y = emmean, x = Duration, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(Type ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit.n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration (ms)", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# n170.all.LinePlot.slide <- {
# n170.all.LinePlot +
# scale_y_reverse(limits = c(0.8, -4.8), breaks = seq(0, -4, -1)) + # set the limit for y axis
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('N170_all.png', n170.all.LinePlot.slide, width = 7, height = 4)
n170.all.LinePlot
#### Contrasts
contra.all.N170 <- contrast(emm.N170.E205, "pairwise", simple = "each", combine = TRUE, adjust = "none")
contra.all.N170
## Category Hemisphere Type Duration contrast estimate SE df t.ratio p.value
## face Left intact . 17 - 200 0.676342 0.2582103 100.06 2.619 0.0102
## house Left intact . 17 - 200 -1.246277 0.2583720 100.31 -4.824 <.0001
## face Right intact . 17 - 200 1.295381 0.2582103 100.06 5.017 <.0001
## house Right intact . 17 - 200 -0.942840 0.2583720 100.31 -3.649 0.0004
## face Left scrambled . 17 - 200 -0.993654 0.2619497 105.98 -3.793 0.0002
## house Left scrambled . 17 - 200 -1.146543 0.2619183 105.93 -4.377 <.0001
## face Right scrambled . 17 - 200 -1.094031 0.2619497 105.98 -4.176 0.0001
## house Right scrambled . 17 - 200 -1.030586 0.2619183 105.93 -3.935 0.0001
## . Left intact 17 face - house -0.714902 0.1954714 116.95 -3.657 0.0004
## . Left intact 200 face - house -2.637521 0.2219773 194.49 -11.882 <.0001
## . Right intact 17 face - house -1.041124 0.1954714 116.95 -5.326 <.0001
## . Right intact 200 face - house -3.279344 0.2219773 194.49 -14.773 <.0001
## . Left scrambled 17 face - house 0.035906 0.2046084 140.40 0.175 0.8609
## . Left scrambled 200 face - house -0.116982 0.2222826 195.56 -0.526 0.5993
## . Right scrambled 17 face - house 0.032712 0.2046084 140.40 0.160 0.8732
## . Right scrambled 200 face - house 0.096157 0.2222826 195.56 0.433 0.6658
## face . intact 17 Left - Right 0.124897 0.2527268 62.97 0.494 0.6229
## face . intact 200 Left - Right 0.743936 0.2735580 86.44 2.719 0.0079
## house . intact 17 Left - Right -0.201325 0.2526993 62.95 -0.797 0.4286
## house . intact 200 Left - Right 0.102112 0.2738871 86.86 0.373 0.7102
## face . scrambled 17 Left - Right -0.198354 0.2598409 70.37 -0.763 0.4478
## face . scrambled 200 Left - Right -0.298730 0.2739983 87.00 -1.090 0.2786
## house . scrambled 17 Left - Right -0.201548 0.2598447 70.37 -0.776 0.4406
## house . scrambled 200 Left - Right -0.085591 0.2739377 86.92 -0.312 0.7554
## face Left . 17 intact - scrambled -0.822638 0.2246563 126.17 -3.662 0.0004
## face Left . 200 intact - scrambled -2.492634 0.2444100 176.75 -10.199 <.0001
## house Left . 17 intact - scrambled -0.071830 0.2246428 126.14 -0.320 0.7497
## house Left . 200 intact - scrambled 0.027904 0.2445602 177.18 0.114 0.9093
## face Right . 17 intact - scrambled -1.145889 0.2246563 126.17 -5.101 <.0001
## face Right . 200 intact - scrambled -3.535300 0.2444100 176.75 -14.465 <.0001
## house Right . 17 intact - scrambled -0.072053 0.2246428 126.14 -0.321 0.7489
## house Right . 200 intact - scrambled -0.159799 0.2445602 177.18 -0.653 0.5143
##
## Degrees-of-freedom method: satterthwaite
confint(contra.all.N170)
contra.all.N170.gen <- contrast(emmeans(lmm.opt.N170.E205, ~ Duration | Hemisphere + Category + Type),
"pairwise")
contra.all.N170.gen[1:2]
## contrast Hemisphere Category Type estimate SE df t.ratio p.value
## 17 - 200 Left face intact 0.6763417 0.2582103 100.06 2.619 0.0102
## 17 - 200 Right face intact 1.2953806 0.2582103 100.06 5.017 <.0001
##
## Degrees-of-freedom method: satterthwaite
confint(contra.all.N170.gen[1:2])
contra.all.N170.gen.hemi <- contrast(emmeans(lmm.opt.N170.E205, ~ Duration + Hemisphere |Category + Type),
interaction = "pairwise")
contra.all.N170.gen.hemi[1]
## Duration_pairwise Hemisphere_pairwise Category Type estimate SE df t.ratio p.value
## 17 - 200 Left - Right face intact -0.6190389 0.2636081 45.58 -2.348 0.0233
##
## Degrees-of-freedom method: satterthwaite
confint(contra.all.N170.gen.hemi[1])
contra.all.N170.spec <- contrast(emmeans(lmm.opt.N170.E205, ~ Category + Duration | Hemisphere + Type),
interaction = "pairwise", adjust = "Bonferroni")
contra.all.N170.spec[c(1,2)]
## Category_pairwise Duration_pairwise Hemisphere Type estimate SE df t.ratio p.value
## face - house 17 - 200 Left intact 1.922618 0.2151955 104.11 8.934 <.0001
## face - house 17 - 200 Right intact 2.238221 0.2151955 104.11 10.401 <.0001
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 2 tests
confint(contra.all.N170.spec[c(1,2)])
contra.all.N170.spec.hemi <- contrast(emmeans(lmm.opt.N170.E205, ~ Category + Duration + Hemisphere | Type),
interaction = "pairwise")
contra.all.N170.spec.hemi
## Type = intact:
## Category_pairwise Duration_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## face - house 17 - 200 Left - Right -0.3156022 0.1910371 70540.74 -1.652 0.0985
##
## Type = scrambled:
## Category_pairwise Duration_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## face - house 17 - 200 Left - Right 0.2163333 0.2099393 70546.20 1.030 0.3028
##
## Degrees-of-freedom method: satterthwaite
confint(contra.all.N170.spec[c(1,2)])
contra.all.N170.gen[1:8]
## contrast Hemisphere Category Type estimate SE df t.ratio p.value
## 17 - 200 Left face intact 0.6763417 0.2582103 100.06 2.619 0.1614
## 17 - 200 Right face intact 1.2953806 0.2582103 100.06 5.017 0.0001
## 17 - 200 Left house intact -1.2462767 0.2583720 100.31 -4.824 0.0001
## 17 - 200 Right house intact -0.9428400 0.2583720 100.31 -3.649 0.0096
## 17 - 200 Left face scrambled -0.9936544 0.2619497 105.98 -3.793 0.0059
## 17 - 200 Right face scrambled -1.0940310 0.2619497 105.98 -4.176 0.0015
## 17 - 200 Left house scrambled -1.1465426 0.2619183 105.93 -4.377 0.0007
## 17 - 200 Right house scrambled -1.0305859 0.2619183 105.93 -3.935 0.0036
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 8 estimates
confint(contra.all.N170.gen[1:8])
contra.all.N170.gen.hemi[1:4]
## Duration_pairwise Hemisphere_pairwise Category Type estimate SE df t.ratio p.value
## 17 - 200 Left - Right face intact -0.6190389 0.2636081 45.58 -2.348 0.0233
## 17 - 200 Left - Right house intact -0.3034367 0.2639238 45.80 -1.150 0.2562
## 17 - 200 Left - Right face scrambled 0.1003766 0.2708823 50.83 0.371 0.7125
## 17 - 200 Left - Right house scrambled -0.1159567 0.2708231 50.78 -0.428 0.6703
##
## Degrees-of-freedom method: satterthwaite
confint(contra.all.N170.gen.hemi[1:4])
contra.all.N170.spec
## Hemisphere = Left, Type = intact:
## Category_pairwise Duration_pairwise estimate SE df t.ratio p.value
## face - house 17 - 200 1.9226184 0.2151955 104.11 8.934 <.0001
##
## Hemisphere = Right, Type = intact:
## Category_pairwise Duration_pairwise estimate SE df t.ratio p.value
## face - house 17 - 200 2.2382207 0.2151955 104.11 10.401 <.0001
##
## Hemisphere = Left, Type = scrambled:
## Category_pairwise Duration_pairwise estimate SE df t.ratio p.value
## face - house 17 - 200 0.1528882 0.2238302 121.85 0.683 0.4959
##
## Hemisphere = Right, Type = scrambled:
## Category_pairwise Duration_pairwise estimate SE df t.ratio p.value
## face - house 17 - 200 -0.0634451 0.2238302 121.85 -0.283 0.7773
##
## Degrees-of-freedom method: satterthwaite
# contra.all.N170.spec.hemi
df.erp.SC.E2 <- {
df.erp.E2 %>%
filter(Type == "intact") %>% # only keep the intact trials
mutate(Confidence = if_else(Response %in% c("11", "55"), "high",
if_else(Response %in% c("12", "54"), "low",
if_else(substring(Response, 2) == "3", "guess", "NA"))),
DuraConf = paste(Duration, Confidence, sep = "_")) %>%
filter(Confidence != "NA",
!(DuraConf %in% c("200_guess", "200_low"))) %>%
mutate(Confidence = as.factor(Confidence),
DuraConf = as.factor(DuraConf))
}
# successive difference coding
df.erp.SC.E2 <- sdif_coding_E205_erp(df.erp.SC.E2)
# only keep the data for amplitudes of the P1 (correct and incorrect erps)
df.erp.P1.SC.E2 <- {
df.erp.SC.E2 %>%
filter(Component == "P1")
}
if (saveCSV) {
output.erp.P1.SC.E2 <- file.path(folder.nesi, "E205_erp_P1_SC.RData")
save(df.erp.P1.SC.E2, file = output.erp.P1.SC.E2)
}
# the structure of this dataset
head(df.erp.P1.SC.E2, 10)
lmm.max.P1.SC.E2 <- lmer(MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 + Hemi_C | Stimuli),
data = df.erp.P1.SC.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa")) # nloptwrap Nelder_Mead
## start par. = 1 0 1 1 0 0 0 1 0 0 1 0 1 fn = 221091.6
## At return
## eval: 1276 fn: 220102.61 par: 0.0378072 0.00819981 0.00000 0.338147 -0.0678826 0.0508097 -0.000956547 0.0773424 0.0621324 0.0254029 0.214512 -0.0111323 2.21822e-06
## boundary (singular) fit: see ?isSingular
print(summary(lmm.max.P1.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category * Hemisphere * DuraConf + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 + Hemi_C | Stimuli)
## Data: df.erp.P1.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 220102.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6020 -0.5787 -0.0155 0.5775 8.4272
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.021529 0.14673
## Hemi_C 0.001013 0.03182 1.00
## SubjCode (Intercept) 1.722174 1.31232
## Cate_C 0.159499 0.39937 -0.66
## Hemi_C 0.790084 0.88887 0.22 0.06
## Cate_Hemi 0.011600 0.10770 -0.03 0.71 -0.14
## Residual 15.061410 3.88090
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.248e+00 2.495e-01 2.730e+01 5.002 2.94e-05 ***
## Category2-1 5.822e-02 9.297e-02 3.683e+01 0.626 0.5350
## Hemisphere2-1 5.166e-02 1.733e-01 2.748e+01 0.298 0.7678
## DuraConf2-1 1.051e-01 6.125e-02 2.983e+04 1.716 0.0862 .
## DuraConf3-2 2.647e-01 6.564e-02 1.776e+04 4.033 5.54e-05 ***
## DuraConf4-3 3.971e-02 6.696e-02 2.288e+04 0.593 0.5532
## Category2-1:Hemisphere2-1 -9.480e-02 8.859e-02 2.180e+02 -1.070 0.2857
## Category2-1:DuraConf2-1 -2.347e-01 1.209e-01 1.212e+04 -1.941 0.0523 .
## Category2-1:DuraConf3-2 3.841e-02 1.263e-01 2.598e+03 0.304 0.7612
## Category2-1:DuraConf4-3 5.656e-02 1.312e-01 6.594e+03 0.431 0.6663
## Hemisphere2-1:DuraConf2-1 1.840e-02 1.205e-01 3.585e+04 0.153 0.8786
## Hemisphere2-1:DuraConf3-2 7.275e-02 1.272e-01 2.436e+04 0.572 0.5673
## Hemisphere2-1:DuraConf4-3 5.993e-02 1.314e-01 2.522e+04 0.456 0.6483
## Category2-1:Hemisphere2-1:DuraConf2-1 1.322e-01 2.334e-01 2.924e+04 0.566 0.5713
## Category2-1:Hemisphere2-1:DuraConf3-2 -1.058e-01 2.330e-01 1.268e+04 -0.454 0.6498
## Category2-1:Hemisphere2-1:DuraConf4-3 -1.807e-01 2.513e-01 2.114e+04 -0.719 0.4721
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
lmm.zcp.P1.SC.E2 <- update(lmm.max.P1.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 + Hemi_C || Stimuli))
## start par. = 1 1 1 1 1 1 fn = 221091.6
## At return
## eval: 509 fn: 220118.27 par: 6.22133e-08 0.0375916 1.47841e-06 0.229002 0.102779 0.338213
## boundary (singular) fit: see ?isSingular
print(summary(lmm.zcp.P1.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C +
## Cate_Hemi || SubjCode) + (1 + Hemi_C || Stimuli) + Category:Hemisphere +
## Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.P1.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 220118.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.604 -0.578 -0.015 0.578 8.456
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Hemi_C 5.830e-14 2.415e-07
## Stimuli.1 (Intercept) 2.129e-02 1.459e-01
## SubjCode Cate_Hemi 3.292e-11 5.738e-06
## SubjCode.1 Hemi_C 7.899e-01 8.888e-01
## SubjCode.2 Cate_C 1.591e-01 3.989e-01
## SubjCode.3 (Intercept) 1.723e+00 1.313e+00
## Residual 1.506e+01 3.881e+00
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.249e+00 2.495e-01 2.730e+01 5.007 2.91e-05 ***
## Category2-1 6.073e-02 9.286e-02 3.677e+01 0.654 0.5172
## Hemisphere2-1 5.095e-02 1.732e-01 2.746e+01 0.294 0.7709
## DuraConf2-1 1.039e-01 6.148e-02 3.568e+04 1.690 0.0911 .
## DuraConf3-2 2.724e-01 6.610e-02 2.887e+04 4.121 3.78e-05 ***
## DuraConf4-3 3.293e-02 6.739e-02 3.365e+04 0.489 0.6251
## Category2-1:Hemisphere2-1 -9.713e-02 8.593e-02 3.917e+04 -1.130 0.2584
## Category2-1:DuraConf2-1 -2.413e-01 1.220e-01 2.111e+04 -1.979 0.0479 *
## Category2-1:DuraConf3-2 6.113e-02 1.290e-01 5.815e+03 0.474 0.6357
## Category2-1:DuraConf4-3 4.158e-02 1.330e-01 1.365e+04 0.313 0.7546
## Hemisphere2-1:DuraConf2-1 2.186e-02 1.205e-01 3.789e+04 0.181 0.8561
## Hemisphere2-1:DuraConf3-2 6.451e-02 1.272e-01 2.990e+04 0.507 0.6121
## Hemisphere2-1:DuraConf4-3 6.541e-02 1.315e-01 3.598e+04 0.497 0.6189
## Category2-1:Hemisphere2-1:DuraConf2-1 1.270e-01 2.332e-01 3.936e+04 0.545 0.5860
## Category2-1:Hemisphere2-1:DuraConf3-2 -1.051e-01 2.325e-01 3.908e+04 -0.452 0.6513
## Category2-1:Hemisphere2-1:DuraConf4-3 -1.768e-01 2.512e-01 3.933e+04 -0.704 0.4814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.zcp.P1.SC.E2))
## $Stimuli
## Importance of components:
## [,1] [,2]
## Standard deviation 0.03759 6.221e-08
## Proportion of Variance 1.00000 0.000e+00
## Cumulative Proportion 1.00000 1.000e+00
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.3382 0.2290 0.10278 1.478e-06
## Proportion of Variance 0.6448 0.2956 0.05955 0.000e+00
## Cumulative Proportion 0.6448 0.9405 1.00000 1.000e+00
lmm.rdc.P1.SC.E2 <- update(lmm.zcp.P1.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C || SubjCode) +
(1 | Stimuli))
## start par. = 1 1 1 1 fn = 220659.7
## At return
## eval: 243 fn: 220118.27 par: 0.0375916 0.229002 0.102779 0.338212
print(summary(lmm.rdc.P1.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C ||
## SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.P1.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 220118.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.604 -0.578 -0.015 0.578 8.456
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli (Intercept) 0.02129 0.1459
## SubjCode Hemi_C 0.78991 0.8888
## SubjCode.1 Cate_C 0.15911 0.3989
## SubjCode.2 (Intercept) 1.72297 1.3126
## Residual 15.06265 3.8811
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.249e+00 2.495e-01 2.730e+01 5.007 2.91e-05 ***
## Category2-1 6.073e-02 9.286e-02 3.677e+01 0.654 0.5172
## Hemisphere2-1 5.095e-02 1.732e-01 2.746e+01 0.294 0.7709
## DuraConf2-1 1.039e-01 6.148e-02 3.568e+04 1.690 0.0911 .
## DuraConf3-2 2.724e-01 6.610e-02 2.887e+04 4.121 3.78e-05 ***
## DuraConf4-3 3.293e-02 6.739e-02 3.365e+04 0.489 0.6251
## Category2-1:Hemisphere2-1 -9.713e-02 8.593e-02 3.917e+04 -1.130 0.2584
## Category2-1:DuraConf2-1 -2.413e-01 1.220e-01 2.111e+04 -1.979 0.0479 *
## Category2-1:DuraConf3-2 6.113e-02 1.290e-01 5.815e+03 0.474 0.6357
## Category2-1:DuraConf4-3 4.158e-02 1.330e-01 1.365e+04 0.313 0.7546
## Hemisphere2-1:DuraConf2-1 2.186e-02 1.205e-01 3.789e+04 0.181 0.8561
## Hemisphere2-1:DuraConf3-2 6.451e-02 1.272e-01 2.990e+04 0.507 0.6121
## Hemisphere2-1:DuraConf4-3 6.541e-02 1.315e-01 3.598e+04 0.497 0.6189
## Category2-1:Hemisphere2-1:DuraConf2-1 1.270e-01 2.332e-01 3.936e+04 0.545 0.5860
## Category2-1:Hemisphere2-1:DuraConf3-2 -1.051e-01 2.325e-01 3.908e+04 -0.452 0.6513
## Category2-1:Hemisphere2-1:DuraConf4-3 -1.768e-01 2.512e-01 3.933e+04 -0.704 0.4814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.P1.SC.E2, lmm.zcp.P1.SC.E2, refit = FALSE)
lmm.etd.P1.SC.E2 <- update(lmm.rdc.P1.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C | SubjCode) +
(1 | Stimuli))
## start par. = 1 1 0 0 1 0 1 fn = 220659.7
## At return
## eval: 531 fn: 220104.36 par: 0.0376818 0.338066 -0.0679951 0.0508008 0.0772716 0.0621638 0.214439
print(summary(lmm.etd.P1.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C |
## SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.P1.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 220104.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6075 -0.5789 -0.0158 0.5770 8.4361
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.02139 0.1462
## SubjCode (Intercept) 1.72148 1.3121
## Cate_C 0.15958 0.3995 -0.66
## Hemi_C 0.78972 0.8887 0.22 0.06
## Residual 15.06256 3.8811
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.247e+00 2.494e-01 2.730e+01 5.001 2.95e-05 ***
## Category2-1 5.783e-02 9.295e-02 3.674e+01 0.622 0.5377
## Hemisphere2-1 5.142e-02 1.732e-01 2.746e+01 0.297 0.7688
## DuraConf2-1 1.059e-01 6.127e-02 3.012e+04 1.728 0.0840 .
## DuraConf3-2 2.622e-01 6.568e-02 1.817e+04 3.992 6.58e-05 ***
## DuraConf4-3 4.173e-02 6.698e-02 2.308e+04 0.623 0.5333
## Category2-1:Hemisphere2-1 -9.602e-02 8.592e-02 3.911e+04 -1.118 0.2638
## Category2-1:DuraConf2-1 -2.375e-01 1.210e-01 1.240e+04 -1.962 0.0497 *
## Category2-1:DuraConf3-2 3.344e-02 1.265e-01 2.649e+03 0.264 0.7916
## Category2-1:DuraConf4-3 6.136e-02 1.313e-01 6.668e+03 0.467 0.6402
## Hemisphere2-1:DuraConf2-1 2.104e-02 1.205e-01 3.756e+04 0.175 0.8614
## Hemisphere2-1:DuraConf3-2 6.857e-02 1.271e-01 2.857e+04 0.539 0.5896
## Hemisphere2-1:DuraConf4-3 6.261e-02 1.314e-01 3.537e+04 0.476 0.6338
## Category2-1:Hemisphere2-1:DuraConf2-1 1.247e-01 2.332e-01 3.933e+04 0.535 0.5928
## Category2-1:Hemisphere2-1:DuraConf3-2 -9.808e-02 2.324e-01 3.898e+04 -0.422 0.6731
## Category2-1:Hemisphere2-1:DuraConf4-3 -1.824e-01 2.511e-01 3.930e+04 -0.726 0.4677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.P1.SC.E2, lmm.etd.P1.SC.E2, refit = FALSE)
lmm.opt.P1.SC.E205 <- lmm.etd.P1.SC.E2
# print(summary(lmm.opt.P1.SC.E2), corr = FALSE)
summary(lmm.opt.P1.SC.E205)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C |
## SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.P1.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 220104.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6075 -0.5789 -0.0158 0.5770 8.4361
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.02139 0.1462
## SubjCode (Intercept) 1.72148 1.3121
## Cate_C 0.15958 0.3995 -0.66
## Hemi_C 0.78972 0.8887 0.22 0.06
## Residual 15.06256 3.8811
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.247e+00 2.494e-01 2.730e+01 5.001 2.95e-05 ***
## Category2-1 5.783e-02 9.295e-02 3.674e+01 0.622 0.5377
## Hemisphere2-1 5.142e-02 1.732e-01 2.746e+01 0.297 0.7688
## DuraConf2-1 1.059e-01 6.127e-02 3.012e+04 1.728 0.0840 .
## DuraConf3-2 2.622e-01 6.568e-02 1.817e+04 3.992 6.58e-05 ***
## DuraConf4-3 4.173e-02 6.698e-02 2.308e+04 0.623 0.5333
## Category2-1:Hemisphere2-1 -9.602e-02 8.592e-02 3.911e+04 -1.118 0.2638
## Category2-1:DuraConf2-1 -2.375e-01 1.210e-01 1.240e+04 -1.962 0.0497 *
## Category2-1:DuraConf3-2 3.344e-02 1.265e-01 2.649e+03 0.264 0.7916
## Category2-1:DuraConf4-3 6.136e-02 1.313e-01 6.668e+03 0.467 0.6402
## Hemisphere2-1:DuraConf2-1 2.104e-02 1.205e-01 3.756e+04 0.175 0.8614
## Hemisphere2-1:DuraConf3-2 6.857e-02 1.271e-01 2.857e+04 0.539 0.5896
## Hemisphere2-1:DuraConf4-3 6.261e-02 1.314e-01 3.537e+04 0.476 0.6338
## Category2-1:Hemisphere2-1:DuraConf2-1 1.247e-01 2.332e-01 3.933e+04 0.535 0.5928
## Category2-1:Hemisphere2-1:DuraConf3-2 -9.808e-02 2.324e-01 3.898e+04 -0.422 0.6731
## Category2-1:Hemisphere2-1:DuraConf4-3 -1.824e-01 2.511e-01 3.930e+04 -0.726 0.4677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
confint(lmm.opt.P1.SC.E205, method = "Wald")
## 2.5 % 97.5 %
## .sig01 NA NA
## .sig02 NA NA
## .sig03 NA NA
## .sig04 NA NA
## .sig05 NA NA
## .sig06 NA NA
## .sig07 NA NA
## .sigma NA NA
## (Intercept) 0.75856142 1.7362659230
## Category2-1 -0.12434217 0.2399977383
## Hemisphere2-1 -0.28803158 0.3908789061
## DuraConf2-1 -0.01421216 0.2259624238
## DuraConf3-2 0.13345668 0.3909005188
## DuraConf4-3 -0.08955721 0.1730168953
## Category2-1:Hemisphere2-1 -0.26442374 0.0723833894
## Category2-1:DuraConf2-1 -0.47465631 -0.0003052446
## Category2-1:DuraConf3-2 -0.21453862 0.2814168415
## Category2-1:DuraConf4-3 -0.19592559 0.3186523455
## Hemisphere2-1:DuraConf2-1 -0.21510182 0.2571774938
## Hemisphere2-1:DuraConf3-2 -0.18060402 0.3177447840
## Hemisphere2-1:DuraConf4-3 -0.19499835 0.3202179117
## Category2-1:Hemisphere2-1:DuraConf2-1 -0.33231516 0.5816903695
## Category2-1:Hemisphere2-1:DuraConf3-2 -0.55367829 0.3575264802
## Category2-1:Hemisphere2-1:DuraConf4-3 -0.67462443 0.3098548097
anova(lmm.opt.P1.SC.E205)
qqplot_lmer(lmm.opt.P1.SC.E205)
emm.P1.SC.E205 <- emmeans(lmm.opt.P1.SC.E205, ~ DuraConf + Category | Hemisphere)
emm.P1.SC.E205
## Hemisphere = Left:
## DuraConf Category emmean SE df lower.CL upper.CL
## 17_guess face 0.9079673 0.2947474 36.58 0.3105221 1.505413
## 17_low face 1.1532358 0.2832600 31.22 0.5756860 1.730786
## 17_high face 1.3398907 0.2802008 29.89 0.7675549 1.912227
## 200_high face 1.2740378 0.2854063 32.19 0.6928164 1.855259
## 17_guess house 1.1592962 0.2318241 32.23 0.6872180 1.631374
## 17_low house 1.1047401 0.2309510 31.76 0.6341687 1.575311
## 17_high house 1.3738720 0.2507271 43.94 0.8685443 1.879200
## 200_high house 1.4605749 0.2369209 35.19 0.9796916 1.941458
##
## Hemisphere = Right:
## DuraConf Category emmean SE df lower.CL upper.CL
## 17_guess face 0.9411258 0.3233773 34.61 0.2843680 1.597884
## 17_low face 1.1450884 0.3129860 30.38 0.5062175 1.783959
## 17_high face 1.4493516 0.3102156 29.31 0.8151845 2.083519
## 200_high face 1.5373008 0.3149735 31.16 0.8950445 2.179557
## 17_guess house 1.0975529 0.2701625 30.76 0.5463754 1.648731
## 17_low house 1.1263785 0.2693767 30.41 0.5765456 1.676211
## 17_high house 1.4150429 0.2868544 39.02 0.8348323 1.995253
## 200_high house 1.4731631 0.2744723 32.78 0.9146025 2.031724
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
# emmip(emm.P1.SC.E205, Category ~ DuraConf | Hemisphere, CIs = TRUE)
duraconf.labels = c("17ms\nguess", "17ms\nlow", "17ms\nhigh", "200ms\nhigh")
p1.conf.LinePlot = {
ggplot(data = data.frame(emm.P1.SC.E205), aes(y = emmean, x = DuraConf, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit.p1) + # set the limit for y axis
scale_x_discrete(labels = duraconf.labels) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Subjective Confidence", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# p1.conf.LinePlot.slide <- {
# p1.conf.LinePlot +
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('P1_conf.png', p1.conf.LinePlot.slide, width = 8.5, height = 4)
p1.conf.LinePlot
# contra.SC.P1 <- contrast(emmeans(lmm.opt.P1.SC.E205, ~ DuraConf + Category), "pairwise", simple = "each", combine = TRUE)
#
# contra.SC.P1
# confint(contra.SC.P1)
contra.SC.P1.gen <- contrast(emmeans(lmm.opt.P1.SC.E205, ~ DuraConf | Hemisphere + Category), #
"pairwise", adjust = "Bonferroni")
contra.SC.P1.gen[1:6]
## contrast Hemisphere Category estimate SE df t.ratio p.value
## 17_guess - 17_low Left face -0.2452685 0.1354348 37115.10 -1.811 0.4209
## 17_guess - 17_high Left face -0.4319234 0.1376417 20687.47 -3.138 0.0102
## 17_guess - 200_high Left face -0.3660704 0.1425063 36491.93 -2.569 0.0613
## 17_low - 17_high Left face -0.1866549 0.1094711 17691.77 -1.705 0.5292
## 17_low - 200_high Left face -0.1208019 0.1171315 37690.70 -1.031 1.0000
## 17_high - 200_high Left face 0.0658529 0.1097015 36231.42 0.600 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
confint(contra.SC.P1.gen[1:6])
contra.SC.P1.gen[7:12]
## contrast Hemisphere Category estimate SE df t.ratio p.value
## 17_guess - 17_low Right face -0.2039625 0.1352749 35346.58 -1.508 0.7897
## 17_guess - 17_high Right face -0.5082257 0.1370905 15779.38 -3.707 0.0013
## 17_guess - 200_high Right face -0.5961750 0.1423218 34614.66 -4.189 0.0002
## 17_low - 17_high Right face -0.3042632 0.1089672 13017.90 -2.792 0.0315
## 17_low - 200_high Right face -0.3922125 0.1170154 36493.11 -3.352 0.0048
## 17_high - 200_high Right face -0.0879492 0.1095568 34348.22 -0.803 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
confint(contra.SC.P1.gen[7:12])
# differences between hemispheres
contra.SC.P1.gen.hemi <- contrast(emmeans(lmm.opt.P1.SC.E205, ~ DuraConf + Hemisphere | Category), #
interaction = "pairwise", adjust = "Bonferroni")
contra.SC.P1.gen.hemi
## Category = face:
## DuraConf_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## 17_guess - 17_low Left - Right -0.04130597 0.1894667 39164.20 -0.218 1.0000
## 17_guess - 17_high Left - Right 0.07630237 0.1884018 34438.11 0.405 1.0000
## 17_guess - 200_high Left - Right 0.23010456 0.1994051 39057.60 1.154 1.0000
## 17_low - 17_high Left - Right 0.11760834 0.1490986 33092.64 0.789 1.0000
## 17_low - 200_high Left - Right 0.27141052 0.1643000 39274.39 1.652 0.5914
## 17_high - 200_high Left - Right 0.15380219 0.1535298 39052.83 1.002 1.0000
##
## Category = house:
## DuraConf_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## 17_guess - 17_low Left - Right 0.08338164 0.1425409 37230.92 0.585 1.0000
## 17_guess - 17_high Left - Right 0.10291407 0.1998271 30438.15 0.515 1.0000
## 17_guess - 200_high Left - Right 0.07433145 0.1591297 39289.66 0.467 1.0000
## 17_low - 17_high Left - Right 0.01953243 0.1926455 35036.26 0.101 1.0000
## 17_low - 200_high Left - Right -0.00905019 0.1571648 39426.36 -0.058 1.0000
## 17_high - 200_high Left - Right -0.02858262 0.2061952 36837.59 -0.139 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
contra.SC.P1.gen.nohemi <- contrast(emmeans(lmm.opt.P1.SC.E205, ~ DuraConf | Category), #
"pairwise", adjust = "Bonferroni")
## NOTE: Results may be misleading due to involvement in interactions
contra.SC.P1.gen.nohemi[1:6]
## contrast Category estimate SE df t.ratio p.value
## 17_guess - 17_low face -0.2246155 0.09667744 31481.61 -2.323 0.1210
## 17_guess - 17_high face -0.4700746 0.09997852 9957.84 -4.702 <.0001
## 17_guess - 200_high face -0.4811227 0.10169157 30213.22 -4.731 <.0001
## 17_low - 17_high face -0.2454590 0.07982037 7934.70 -3.075 0.0127
## 17_low - 200_high face -0.2565072 0.08341207 33288.82 -3.075 0.0126
## 17_high - 200_high face -0.0110482 0.07826687 29623.74 -0.141 1.0000
##
## Results are averaged over the levels of: Hemisphere
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
confint(contra.SC.P1.gen.nohemi[1:6])
contra.SC.P1.spec <- contrast(emmeans(lmm.opt.P1.SC.E205, ~ Category + DuraConf | Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni")
contra.SC.P1.spec
## Hemisphere = Left:
## Category_pairwise DuraConf_pairwise estimate SE df t.ratio p.value
## face - house 17_guess - 17_low -0.29982458 0.1679815 24749.21 -1.785 0.4458
## face - house 17_guess - 17_high -0.21734752 0.1920037 5076.78 -1.132 1.0000
## face - house 17_guess - 200_high -0.06479173 0.1800928 31977.14 -0.360 1.0000
## face - house 17_low - 17_high 0.08247706 0.1716760 7492.40 0.480 1.0000
## face - house 17_low - 200_high 0.23503285 0.1617098 38100.50 1.453 0.8767
## face - house 17_high - 200_high 0.15255578 0.1814463 15301.15 0.841 1.0000
##
## Hemisphere = Right:
## Category_pairwise DuraConf_pairwise estimate SE df t.ratio p.value
## face - house 17_guess - 17_low -0.17513697 0.1680863 25502.54 -1.042 1.0000
## face - house 17_guess - 17_high -0.19073582 0.1923894 5541.48 -0.991 1.0000
## face - house 17_guess - 200_high -0.22056484 0.1800024 31001.30 -1.225 1.0000
## face - house 17_low - 17_high -0.01559884 0.1719311 8082.12 -0.091 1.0000
## face - house 17_low - 200_high -0.04542787 0.1616989 37871.30 -0.281 1.0000
## face - house 17_high - 200_high -0.02982903 0.1818779 17940.90 -0.164 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
confint(contra.SC.P1.spec)
contra.SC.P1.spec.hemi <- contrast(emmeans(lmm.opt.P1.SC.E205, ~ Category + DuraConf + Hemisphere ), #
interaction = "pairwise", adjust = "Bonferroni")
contra.SC.P1.spec.hemi
## Category_pairwise DuraConf_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## face - house 17_guess - 17_low Left - Right -0.1246876 0.2331690 39329.59 -0.535 1.0000
## face - house 17_guess - 17_high Left - Right -0.0266117 0.2573193 39170.26 -0.103 1.0000
## face - house 17_guess - 200_high Left - Right 0.1557731 0.2514099 39451.95 0.620 1.0000
## face - house 17_low - 17_high Left - Right 0.0980759 0.2324545 38979.93 0.422 1.0000
## face - house 17_low - 200_high Left - Right 0.2804607 0.2272674 39363.21 1.234 1.0000
## face - house 17_high - 200_high Left - Right 0.1823848 0.2511473 39297.25 0.726 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
# only keep the data for amplitudes of the N170 (all erps)
df.erp.N170.SC.E2 <- {
df.erp.SC.E2 %>%
filter(Component == "N170")
}
if (saveCSV) {
output.erp.N170.SC.E2 <- file.path(folder.nesi, "E205_erp_N170_SC.RData")
save(df.erp.N170.SC.E2, file = output.erp.N170.SC.E2)
}
# the structure of this dataset
head(df.erp.N170.SC.E2, 10)
lmm.max.N170.SC.E2 <- lmer(MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 + Hemi_C | Stimuli),
data = df.erp.N170.SC.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa")) # nloptwrap Nelder_Mead
## start par. = 1 0 1 1 0 0 0 1 0 0 1 0 1 fn = 222383.5
## At return
## eval: 1493 fn: 221500.89 par: 0.0377450 0.00592574 1.37979e-07 0.468918 -0.124988 -0.00676750 -0.0181408 0.172581 -0.0131746 -0.0558577 0.351848 -0.213332 0.184848
## boundary (singular) fit: see ?isSingular
print(summary(lmm.max.N170.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category * Hemisphere * DuraConf + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 + Hemi_C | Stimuli)
## Data: df.erp.N170.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 221500.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3165 -0.5946 0.0000 0.5895 6.7424
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 2.217e-02 0.14888
## Hemi_C 5.463e-04 0.02337 1.00
## SubjCode (Intercept) 3.421e+00 1.84957
## Cate_C 7.064e-01 0.84049 -0.59
## Hemi_C 1.929e+00 1.38904 -0.02 -0.02
## Cate_Hemi 1.293e+00 1.13723 -0.06 -0.12 -0.73
## Residual 1.556e+01 3.94434
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.996e+00 3.506e-01 2.714e+01 -5.693 4.68e-06 ***
## Category2-1 1.217e+00 1.682e-01 2.946e+01 7.236 5.24e-08 ***
## Hemisphere2-1 -1.775e-01 2.661e-01 2.726e+01 -0.667 0.51039
## DuraConf2-1 -4.935e-01 6.263e-02 3.817e+04 -7.880 3.35e-15 ***
## DuraConf3-2 -1.805e-01 6.747e-02 3.568e+04 -2.675 0.00747 **
## DuraConf4-3 3.715e-01 6.865e-02 3.638e+04 5.412 6.29e-08 ***
## Category2-1:Hemisphere2-1 3.189e-01 2.322e-01 2.856e+01 1.373 0.18036
## Category2-1:DuraConf2-1 5.042e-01 1.248e-01 3.307e+04 4.039 5.37e-05 ***
## Category2-1:DuraConf3-2 3.517e-01 1.335e-01 1.716e+04 2.635 0.00843 **
## Category2-1:DuraConf4-3 1.770e+00 1.365e-01 2.594e+04 12.971 < 2e-16 ***
## Hemisphere2-1:DuraConf2-1 -1.560e-03 1.242e-01 3.044e+04 -0.013 0.98998
## Hemisphere2-1:DuraConf3-2 -2.377e-01 1.329e-01 1.699e+04 -1.789 0.07358 .
## Hemisphere2-1:DuraConf4-3 -1.685e-01 1.354e-01 1.950e+04 -1.244 0.21353
## Category2-1:Hemisphere2-1:DuraConf2-1 2.736e-01 2.465e-01 1.737e+04 1.110 0.26707
## Category2-1:Hemisphere2-1:DuraConf3-2 -5.786e-02 2.597e-01 4.298e+03 -0.223 0.82373
## Category2-1:Hemisphere2-1:DuraConf4-3 2.674e-01 2.675e-01 8.591e+03 1.000 0.31755
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
lmm.zcp.N170.SC.E2 <- update(lmm.max.N170.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 + Hemi_C || Stimuli))
## start par. = 1 1 1 1 1 1 fn = 222383.5
## At return
## eval: 564 fn: 221530.50 par: 0.00000 0.0376425 0.288393 0.352222 0.212488 0.469156
## boundary (singular) fit: see ?isSingular
print(summary(lmm.zcp.N170.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C +
## Cate_Hemi || SubjCode) + (1 + Hemi_C || Stimuli) + Category:Hemisphere +
## Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.N170.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 221530.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3076 -0.5935 -0.0001 0.5902 6.7081
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Hemi_C 0.00000 0.0000
## Stimuli.1 (Intercept) 0.02205 0.1485
## SubjCode Cate_Hemi 1.29398 1.1375
## SubjCode.1 Hemi_C 1.93014 1.3893
## SubjCode.2 Cate_C 0.70246 0.8381
## SubjCode.3 (Intercept) 3.42446 1.8505
## Residual 15.55810 3.9444
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.995e+00 3.508e-01 2.714e+01 -5.687 4.76e-06 ***
## Category2-1 1.217e+00 1.677e-01 2.946e+01 7.253 5.00e-08 ***
## Hemisphere2-1 -1.757e-01 2.662e-01 2.726e+01 -0.660 0.51480
## DuraConf2-1 -4.958e-01 6.270e-02 3.886e+04 -7.909 2.67e-15 ***
## DuraConf3-2 -1.785e-01 6.761e-02 3.821e+04 -2.640 0.00829 **
## DuraConf4-3 3.698e-01 6.879e-02 3.855e+04 5.376 7.66e-08 ***
## Category2-1:Hemisphere2-1 3.189e-01 2.323e-01 2.856e+01 1.373 0.18044
## Category2-1:DuraConf2-1 5.018e-01 1.251e-01 3.672e+04 4.011 6.05e-05 ***
## Category2-1:DuraConf3-2 3.702e-01 1.343e-01 2.615e+04 2.757 0.00584 **
## Category2-1:DuraConf4-3 1.759e+00 1.371e-01 3.347e+04 12.835 < 2e-16 ***
## Hemisphere2-1:DuraConf2-1 3.713e-03 1.248e-01 3.817e+04 0.030 0.97626
## Hemisphere2-1:DuraConf3-2 -2.309e-01 1.342e-01 3.377e+04 -1.721 0.08521 .
## Hemisphere2-1:DuraConf4-3 -1.775e-01 1.368e-01 3.704e+04 -1.298 0.19430
## Category2-1:Hemisphere2-1:DuraConf2-1 2.960e-01 2.486e-01 3.101e+04 1.191 0.23377
## Category2-1:Hemisphere2-1:DuraConf3-2 -3.908e-02 2.652e-01 1.341e+04 -0.147 0.88287
## Category2-1:Hemisphere2-1:DuraConf4-3 2.456e-01 2.718e-01 2.410e+04 0.904 0.36614
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(lmm.max.N170.SC.E2, lmm.zcp.N170.SC.E2, refit = FALSE)
summary(rePCA(lmm.zcp.N170.SC.E2))
## $Stimuli
## Importance of components:
## [,1] [,2]
## Standard deviation 0.03764 0
## Proportion of Variance 1.00000 0
## Cumulative Proportion 1.00000 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.4692 0.3522 0.2884 0.21249
## Proportion of Variance 0.4658 0.2626 0.1760 0.09556
## Cumulative Proportion 0.4658 0.7284 0.9044 1.00000
lmm.rdc.N170.SC.E2 <- update(lmm.zcp.N170.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 | Stimuli))
## start par. = 1 1 1 1 1 fn = 222042.5
## At return
## eval: 456 fn: 221530.50 par: 0.0376425 0.288392 0.352219 0.212488 0.469155
print(summary(lmm.rdc.N170.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C +
## Cate_Hemi || SubjCode) + (1 | Stimuli) + Category:Hemisphere +
## Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.N170.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 221530.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3076 -0.5935 -0.0001 0.5902 6.7081
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli (Intercept) 0.02205 0.1485
## SubjCode Cate_Hemi 1.29397 1.1375
## SubjCode.1 Hemi_C 1.93011 1.3893
## SubjCode.2 Cate_C 0.70247 0.8381
## SubjCode.3 (Intercept) 3.42443 1.8505
## Residual 15.55810 3.9444
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.995e+00 3.508e-01 2.714e+01 -5.687 4.75e-06 ***
## Category2-1 1.217e+00 1.677e-01 2.946e+01 7.253 5.00e-08 ***
## Hemisphere2-1 -1.757e-01 2.662e-01 2.726e+01 -0.660 0.51479
## DuraConf2-1 -4.958e-01 6.270e-02 3.886e+04 -7.909 2.67e-15 ***
## DuraConf3-2 -1.785e-01 6.761e-02 3.821e+04 -2.640 0.00829 **
## DuraConf4-3 3.698e-01 6.879e-02 3.855e+04 5.376 7.66e-08 ***
## Category2-1:Hemisphere2-1 3.189e-01 2.323e-01 2.856e+01 1.373 0.18044
## Category2-1:DuraConf2-1 5.018e-01 1.251e-01 3.672e+04 4.011 6.05e-05 ***
## Category2-1:DuraConf3-2 3.702e-01 1.343e-01 2.615e+04 2.757 0.00584 **
## Category2-1:DuraConf4-3 1.759e+00 1.371e-01 3.347e+04 12.835 < 2e-16 ***
## Hemisphere2-1:DuraConf2-1 3.713e-03 1.248e-01 3.817e+04 0.030 0.97626
## Hemisphere2-1:DuraConf3-2 -2.309e-01 1.342e-01 3.377e+04 -1.721 0.08521 .
## Hemisphere2-1:DuraConf4-3 -1.775e-01 1.368e-01 3.704e+04 -1.298 0.19430
## Category2-1:Hemisphere2-1:DuraConf2-1 2.960e-01 2.486e-01 3.101e+04 1.191 0.23377
## Category2-1:Hemisphere2-1:DuraConf3-2 -3.908e-02 2.652e-01 1.341e+04 -0.147 0.88287
## Category2-1:Hemisphere2-1:DuraConf4-3 2.456e-01 2.718e-01 2.410e+04 0.904 0.36614
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.N170.SC.E2, lmm.zcp.N170.SC.E2, refit = FALSE)
lmm.etd.N170.SC.E2 <- update(lmm.rdc.N170.SC.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraConf +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 | Stimuli))
## start par. = 1 1 0 0 0 1 0 0 1 0 1 fn = 222042.5
## At return
## eval: 1134 fn: 221501.14 par: 0.0376527 0.468917 -0.124981 -0.00680051 -0.0181522 0.172586 -0.0130994 -0.0558253 0.351828 -0.213356 0.184805
print(summary(lmm.rdc.N170.SC.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C +
## Cate_Hemi || SubjCode) + (1 | Stimuli) + Category:Hemisphere +
## Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.N170.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 221530.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3076 -0.5935 -0.0001 0.5902 6.7081
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli (Intercept) 0.02205 0.1485
## SubjCode Cate_Hemi 1.29397 1.1375
## SubjCode.1 Hemi_C 1.93011 1.3893
## SubjCode.2 Cate_C 0.70247 0.8381
## SubjCode.3 (Intercept) 3.42443 1.8505
## Residual 15.55810 3.9444
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.995e+00 3.508e-01 2.714e+01 -5.687 4.75e-06 ***
## Category2-1 1.217e+00 1.677e-01 2.946e+01 7.253 5.00e-08 ***
## Hemisphere2-1 -1.757e-01 2.662e-01 2.726e+01 -0.660 0.51479
## DuraConf2-1 -4.958e-01 6.270e-02 3.886e+04 -7.909 2.67e-15 ***
## DuraConf3-2 -1.785e-01 6.761e-02 3.821e+04 -2.640 0.00829 **
## DuraConf4-3 3.698e-01 6.879e-02 3.855e+04 5.376 7.66e-08 ***
## Category2-1:Hemisphere2-1 3.189e-01 2.323e-01 2.856e+01 1.373 0.18044
## Category2-1:DuraConf2-1 5.018e-01 1.251e-01 3.672e+04 4.011 6.05e-05 ***
## Category2-1:DuraConf3-2 3.702e-01 1.343e-01 2.615e+04 2.757 0.00584 **
## Category2-1:DuraConf4-3 1.759e+00 1.371e-01 3.347e+04 12.835 < 2e-16 ***
## Hemisphere2-1:DuraConf2-1 3.713e-03 1.248e-01 3.817e+04 0.030 0.97626
## Hemisphere2-1:DuraConf3-2 -2.309e-01 1.342e-01 3.377e+04 -1.721 0.08521 .
## Hemisphere2-1:DuraConf4-3 -1.775e-01 1.368e-01 3.704e+04 -1.298 0.19430
## Category2-1:Hemisphere2-1:DuraConf2-1 2.960e-01 2.486e-01 3.101e+04 1.191 0.23377
## Category2-1:Hemisphere2-1:DuraConf3-2 -3.908e-02 2.652e-01 1.341e+04 -0.147 0.88287
## Category2-1:Hemisphere2-1:DuraConf4-3 2.456e-01 2.718e-01 2.410e+04 0.904 0.36614
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.N170.SC.E2, lmm.etd.N170.SC.E2, refit = FALSE)
lmm.opt.N170.SC.E205 <- lmm.etd.N170.SC.E2
print(summary(lmm.opt.N170.SC.E205), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraConf + (1 + Cate_C + Hemi_C +
## Cate_Hemi | SubjCode) + (1 | Stimuli) + Category:Hemisphere +
## Category:DuraConf + Hemisphere:DuraConf + Category:Hemisphere:DuraConf
## Data: df.erp.N170.SC.E2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 221501.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3165 -0.5946 -0.0004 0.5897 6.7418
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.02206 0.1485
## SubjCode (Intercept) 3.42096 1.8496
## Cate_C 0.70643 0.8405 -0.59
## Hemi_C 1.92921 1.3890 -0.02 -0.02
## Cate_Hemi 1.29318 1.1372 -0.06 -0.12 -0.73
## Residual 15.55805 3.9444
## Number of obs: 39600, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.996e+00 3.506e-01 2.714e+01 -5.693 4.68e-06 ***
## Category2-1 1.217e+00 1.681e-01 2.945e+01 7.236 5.24e-08 ***
## Hemisphere2-1 -1.775e-01 2.661e-01 2.725e+01 -0.667 0.51031
## DuraConf2-1 -4.936e-01 6.263e-02 3.817e+04 -7.882 3.31e-15 ***
## DuraConf3-2 -1.805e-01 6.747e-02 3.568e+04 -2.674 0.00749 **
## DuraConf4-3 3.715e-01 6.865e-02 3.638e+04 5.411 6.29e-08 ***
## Category2-1:Hemisphere2-1 3.184e-01 2.321e-01 2.853e+01 1.372 0.18085
## Category2-1:DuraConf2-1 5.042e-01 1.248e-01 3.307e+04 4.040 5.36e-05 ***
## Category2-1:DuraConf3-2 3.517e-01 1.335e-01 1.717e+04 2.634 0.00844 **
## Category2-1:DuraConf4-3 1.770e+00 1.365e-01 2.596e+04 12.970 < 2e-16 ***
## Hemisphere2-1:DuraConf2-1 -3.097e-03 1.242e-01 3.044e+04 -0.025 0.98010
## Hemisphere2-1:DuraConf3-2 -2.384e-01 1.329e-01 1.699e+04 -1.794 0.07280 .
## Hemisphere2-1:DuraConf4-3 -1.673e-01 1.354e-01 1.950e+04 -1.235 0.21683
## Category2-1:Hemisphere2-1:DuraConf2-1 2.733e-01 2.465e-01 1.737e+04 1.109 0.26762
## Category2-1:Hemisphere2-1:DuraConf3-2 -5.871e-02 2.597e-01 4.297e+03 -0.226 0.82118
## Category2-1:Hemisphere2-1:DuraConf4-3 2.686e-01 2.675e-01 8.591e+03 1.004 0.31539
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.opt.N170.SC.E205)
qqplot_lmer(lmm.opt.N170.SC.E205)
emm.N170.SC.E205 <- emmeans(lmm.opt.N170.SC.E205, ~ Category + DuraConf | Hemisphere)
emm.N170.SC.E205
## Hemisphere = Left:
## Category DuraConf emmean SE df lower.CL upper.CL
## face 17_guess -1.712560 0.4538844 30.69 -2.638645 -0.7864760
## house 17_guess -1.530219 0.3427950 29.18 -2.231130 -0.8293078
## face 17_low -2.388422 0.4459205 28.59 -3.300995 -1.4758496
## house 17_low -1.838474 0.3420131 28.91 -2.538059 -1.1388892
## face 17_high -2.640200 0.4439514 28.09 -3.549459 -1.7309409
## house 17_high -1.709212 0.3580465 34.68 -2.436328 -0.9820953
## face 200_high -3.003007 0.4469683 28.87 -3.917342 -2.0886717
## house 200_high -0.436159 0.3459560 30.28 -1.142425 0.2701067
##
## Hemisphere = Right:
## Category DuraConf emmean SE df lower.CL upper.CL
## face 17_guess -1.764593 0.4537364 30.67 -2.690398 -0.8387881
## house 17_guess -1.506561 0.3260744 29.31 -2.173154 -0.8399688
## face 17_low -2.580185 0.4458010 28.58 -3.492528 -1.6678427
## house 17_low -1.681280 0.3253093 29.04 -2.346573 -1.0159872
## face 17_high -3.041003 0.4438271 28.08 -3.950024 -2.1319814
## house 17_high -1.819765 0.3415646 35.22 -2.513020 -1.1265102
## face 200_high -3.705354 0.4468824 28.87 -4.619515 -2.7911915
## house 200_high -0.579688 0.3295171 30.58 -1.252118 0.0927434
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
# emmip(emm.N170.SC.E205, Category ~ DuraConf | Hemisphere, CIs = TRUE) +
# scale_y_continuous(trans = "reverse")
duraconf.labels = c("17ms\nguess", "17ms\nlow", "17ms\nhigh", "200ms\nhigh")
n170.conf.LinePlot = {
ggplot(data = data.frame(emm.N170.SC.E205), aes(y = emmean, x = DuraConf, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit.n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
scale_x_discrete(labels = duraconf.labels) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Subjective Confidence", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# n170.conf.LinePlot.slide <- {
# n170.conf.LinePlot +
# scale_y_reverse(limits = c(0.8, -4.8), breaks = seq(0, -4, -1)) + # set the limit for y axis
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('N170_confi.png', n170.conf.LinePlot.slide, width = 8.5, height = 4)
n170.conf.LinePlot
# contra.N170.resp <- contrast(emm.N170.SC.E205, "pairwise", simple = "each", combine = FALSE, adjust = "Bonferroni")
#
# contra.N170.resp
# confint(contra.N170.resp)
contra.SC.N170.gen <- contrast(emmeans(lmm.opt.N170.SC.E205, ~ DuraConf | Hemisphere + Category), #
"pairwise", adjust = "Bonferroni")
contra.SC.N170.gen[1:6]
## contrast Hemisphere Category estimate SE df t.ratio p.value
## 17_guess - 17_low Left face 0.6758617 0.1395425 38713.17 4.843 <.0001
## 17_guess - 17_high Left face 0.9276398 0.1460077 29514.74 6.353 <.0001
## 17_guess - 200_high Left face 1.2904465 0.1469225 38302.56 8.783 <.0001
## 17_low - 17_high Left face 0.2517781 0.1169049 27206.82 2.154 0.1876
## 17_low - 200_high Left face 0.6145848 0.1203932 38772.86 5.105 <.0001
## 17_high - 200_high Left face 0.3628068 0.1131308 38093.93 3.207 0.0081
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
confint(contra.SC.N170.gen[1:6])
contra.SC.N170.gen[7:12]
## contrast Hemisphere Category estimate SE df t.ratio p.value
## 17_guess - 17_low Right face 0.8155921 0.1393891 37587.05 5.851 <.0001
## 17_guess - 17_high Right face 1.2764096 0.1454448 22943.16 8.776 <.0001
## 17_guess - 200_high Right face 1.9407603 0.1467337 36887.27 13.226 <.0001
## 17_low - 17_high Right face 0.4608175 0.1163862 20095.69 3.959 0.0005
## 17_low - 200_high Right face 1.1251682 0.1202709 37878.28 9.355 <.0001
## 17_high - 200_high Right face 0.6643506 0.1129761 36540.10 5.880 <.0001
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
confint(contra.SC.N170.gen[7:12])
contra.SC.N170.gen.hemi <- contrast(emmeans(lmm.opt.N170.SC.E205, ~ DuraConf + Hemisphere | Category), #
interaction = "pairwise", adjust = "Bonferroni")
contra.SC.N170.gen.hemi[1:6]
## DuraConf_pairwise Hemisphere_pairwise Category estimate SE df t.ratio p.value
## 17_guess - 17_low Left - Right face -0.1397304 0.1969419 37121.28 -0.710 1.0000
## 17_guess - 17_high Left - Right face -0.3487699 0.2052731 21129.08 -1.699 0.5359
## 17_guess - 200_high Left - Right face -0.6503138 0.2073552 36333.86 -3.136 0.0103
## 17_low - 17_high Left - Right face -0.2090395 0.1642874 18263.35 -1.272 1.0000
## 17_low - 200_high Left - Right face -0.5105834 0.1700071 37522.47 -3.003 0.0160
## 17_high - 200_high Left - Right face -0.3015439 0.1596700 35954.30 -1.889 0.3538
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
confint(contra.SC.N170.gen.hemi[1:6])
contra.SC.N170.spec <- contrast(emmeans(lmm.opt.N170.SC.E205, ~ Category + DuraConf | Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni")
contra.SC.N170.spec[1:6]
## Category_pairwise DuraConf_pairwise Hemisphere estimate SE df t.ratio p.value
## face - house 17_guess - 17_low Left 0.3676065 0.1758644 29293.02 2.090 0.2196
## face - house 17_guess - 17_high Left 0.7486470 0.2128595 8673.07 3.517 0.0026
## face - house 17_guess - 200_high Left 2.3845062 0.1868966 34616.75 12.758 <.0001
## face - house 17_low - 17_high Left 0.3810405 0.1874598 11931.61 2.033 0.2527
## face - house 17_low - 200_high Left 2.0168996 0.1659301 38713.15 12.155 <.0001
## face - house 17_high - 200_high Left 1.6358591 0.1919170 19722.87 8.524 <.0001
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
confint(contra.SC.N170.spec[1:6])
contra.SC.N170.spec[7:12]
## Category_pairwise DuraConf_pairwise Hemisphere estimate SE df t.ratio p.value
## face - house 17_guess - 17_low Right 0.6408733 0.1749512 20554.15 3.663 0.0015
## face - house 17_guess - 17_high Right 0.9632057 0.2094177 3906.36 4.599 <.0001
## face - house 17_guess - 200_high Right 2.8676340 0.1862807 29077.10 15.394 <.0001
## face - house 17_low - 17_high Right 0.3223323 0.1850158 5675.07 1.742 0.4892
## face - house 17_low - 200_high Right 2.2267607 0.1657006 37310.22 13.438 <.0001
## face - house 17_high - 200_high Right 1.9044283 0.1902512 11186.98 10.010 <.0001
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
confint(contra.SC.N170.spec[7:12])
contra.SC.N170.spec.hemi <- contrast(emmeans(lmm.opt.N170.SC.E205, ~ Category + DuraConf + Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni")
contra.SC.N170.spec.hemi
## Category_pairwise DuraConf_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## face - house 17_guess - 17_low Left - Right -0.2732668 0.2464975 17365.60 -1.109 1.0000
## face - house 17_guess - 17_high Left - Right -0.2145586 0.2931619 2919.40 -0.732 1.0000
## face - house 17_guess - 200_high Left - Right -0.4831278 0.2629621 27102.36 -1.837 0.3971
## face - house 17_low - 17_high Left - Right 0.0587082 0.2597280 4296.86 0.226 1.0000
## face - house 17_low - 200_high Left - Right -0.2098610 0.2341167 36563.69 -0.896 1.0000
## face - house 17_high - 200_high Left - Right -0.2685692 0.2674890 8591.17 -1.004 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 6 tests
tmpBlock.count <- {
clean.beha.205 %>%
filter(Type == "intact") %>%
group_by(SubjCode, Category, Duration, blockID) %>%
summarize(Count_sum = n())
}
avg.resp.long.E205.Block <- {
clean.beha.205 %>%
filter(Type == "intact") %>%
group_by(SubjCode, Category, Duration, blockID, Resp) %>%
summarize(Count = n()) %>%
right_join(tmpBlock.count, by = c("SubjCode", "blockID", "Category", "Duration")) %>%
mutate(RespRate = Count / Count_sum) %>%
ungroup()
}
# descriptive statistics of accuracy for plotting
desc.resp.E205.Block <- {
avg.resp.long.E205.Block %>%
group_by(blockID, Category, Duration, Resp) %>%
summarize(Mean = mean(RespRate)
, N = n()
, SD = if_else(N > 1, sd(RespRate), 0)
, SE = SD/sqrt(N)
, SE.lo = Mean - SE, SE.hi = Mean + SE
, CI = SE * qt(0.975,N)
, Median = median(RespRate)
, Lower = Mean - SD, Upper = Mean + SD # for rainclound plot
) %>%
ungroup() %>%
add_row(blockID = 2, Category = "house", Duration = 200, Resp = 2, Mean = 0, N = 0, SD = 0, SE = 0, SE.lo = 0,
SE.hi = 0, CI = 0, Median = 0, Lower = 0, Upper = 0)
}
avg.resp.wide.E205.Block <- {
avg.resp.long.E205.Block %>%
mutate(Conditions = paste(blockID, Category, Duration, Resp, sep = "_")) %>%
select(SubjCode, Conditions, RespRate) %>%
spread(Conditions, RespRate)
}
block.high <- {
avg.resp.long.E205.Block %>%
filter(Category == "face", Duration == 17, Resp == 1, SubjCode != "P513")
}
t.test(filter(block.high, blockID == 1)$RespRate, filter(block.high, blockID == 2)$RespRate, paired = T)
##
## Paired t-test
##
## data: filter(block.high, blockID == 1)$RespRate and filter(block.high, blockID == 2)$RespRate
## t = -0.64323, df = 26, p-value = 0.5257
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.10424318 0.05455182
## sample estimates:
## mean of the differences
## -0.02484568
knitr::kable(desc.resp.E205.Block, digits = 4)
| blockID | Category | Duration | Resp | Mean | N | SD | SE | SE.lo | SE.hi | CI | Median | Lower | Upper |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | face | 17 | 1 | 0.4628 | 28 | 0.3419 | 0.0646 | 0.3982 | 0.5274 | 0.1324 | 0.3812 | 0.1209 | 0.8047 |
| 1 | face | 17 | 2 | 0.3192 | 28 | 0.2240 | 0.0423 | 0.2769 | 0.3615 | 0.0867 | 0.3354 | 0.0952 | 0.5432 |
| 1 | face | 17 | 3 | 0.1737 | 28 | 0.1856 | 0.0351 | 0.1386 | 0.2087 | 0.0719 | 0.0958 | -0.0120 | 0.3593 |
| 1 | face | 17 | 4 | 0.0472 | 25 | 0.0638 | 0.0128 | 0.0344 | 0.0599 | 0.0263 | 0.0208 | -0.0166 | 0.1109 |
| 1 | face | 17 | 5 | 0.0156 | 4 | 0.0202 | 0.0101 | 0.0055 | 0.0257 | 0.0281 | 0.0062 | -0.0046 | 0.0359 |
| 1 | house | 17 | 1 | 0.0120 | 17 | 0.0107 | 0.0026 | 0.0094 | 0.0146 | 0.0055 | 0.0083 | 0.0013 | 0.0227 |
| 1 | house | 17 | 2 | 0.0658 | 28 | 0.0806 | 0.0152 | 0.0505 | 0.0810 | 0.0312 | 0.0417 | -0.0148 | 0.1464 |
| 1 | house | 17 | 3 | 0.4188 | 28 | 0.2770 | 0.0524 | 0.3664 | 0.4711 | 0.1072 | 0.3542 | 0.1417 | 0.6958 |
| 1 | house | 17 | 4 | 0.3890 | 27 | 0.2334 | 0.0449 | 0.3441 | 0.4340 | 0.0922 | 0.4167 | 0.1556 | 0.6225 |
| 1 | house | 17 | 5 | 0.1961 | 19 | 0.2386 | 0.0547 | 0.1413 | 0.2508 | 0.1146 | 0.1042 | -0.0425 | 0.4347 |
| 2 | face | 17 | 1 | 0.5046 | 27 | 0.3326 | 0.0640 | 0.4406 | 0.5686 | 0.1313 | 0.4750 | 0.1721 | 0.8372 |
| 2 | face | 17 | 2 | 0.3786 | 24 | 0.2660 | 0.0543 | 0.3243 | 0.4329 | 0.1121 | 0.3875 | 0.1126 | 0.6446 |
| 2 | face | 17 | 3 | 0.1300 | 25 | 0.1197 | 0.0239 | 0.1061 | 0.1539 | 0.0493 | 0.0750 | 0.0103 | 0.2497 |
| 2 | face | 17 | 4 | 0.0863 | 21 | 0.0828 | 0.0181 | 0.0682 | 0.1044 | 0.0376 | 0.0500 | 0.0035 | 0.1691 |
| 2 | face | 17 | 5 | 0.0225 | 10 | 0.0184 | 0.0058 | 0.0167 | 0.0283 | 0.0130 | 0.0125 | 0.0041 | 0.0409 |
| 2 | face | 200 | 1 | 0.9768 | 28 | 0.0320 | 0.0061 | 0.9707 | 0.9828 | 0.0124 | 0.9875 | 0.9447 | 1.0088 |
| 2 | face | 200 | 2 | 0.0400 | 10 | 0.0403 | 0.0127 | 0.0273 | 0.0527 | 0.0284 | 0.0188 | -0.0003 | 0.0803 |
| 2 | face | 200 | 3 | 0.0208 | 9 | 0.0125 | 0.0042 | 0.0167 | 0.0250 | 0.0094 | 0.0125 | 0.0083 | 0.0333 |
| 2 | face | 200 | 4 | 0.0125 | 1 | 0.0000 | 0.0000 | 0.0125 | 0.0125 | 0.0000 | 0.0125 | 0.0125 | 0.0125 |
| 2 | face | 200 | 5 | 0.0125 | 4 | 0.0000 | 0.0000 | 0.0125 | 0.0125 | 0.0000 | 0.0125 | 0.0125 | 0.0125 |
| 2 | house | 17 | 1 | 0.0156 | 4 | 0.0062 | 0.0031 | 0.0125 | 0.0188 | 0.0087 | 0.0125 | 0.0094 | 0.0219 |
| 2 | house | 17 | 2 | 0.0354 | 18 | 0.0412 | 0.0097 | 0.0257 | 0.0451 | 0.0204 | 0.0250 | -0.0058 | 0.0766 |
| 2 | house | 17 | 3 | 0.3032 | 27 | 0.2482 | 0.0478 | 0.2555 | 0.3510 | 0.0980 | 0.3000 | 0.0550 | 0.5515 |
| 2 | house | 17 | 4 | 0.4696 | 28 | 0.2377 | 0.0449 | 0.4247 | 0.5146 | 0.0920 | 0.4562 | 0.2319 | 0.7074 |
| 2 | house | 17 | 5 | 0.2208 | 27 | 0.2831 | 0.0545 | 0.1664 | 0.2753 | 0.1118 | 0.0750 | -0.0622 | 0.5039 |
| 2 | house | 200 | 1 | 0.0188 | 6 | 0.0105 | 0.0043 | 0.0145 | 0.0230 | 0.0104 | 0.0125 | 0.0083 | 0.0292 |
| 2 | house | 200 | 3 | 0.0150 | 10 | 0.0079 | 0.0025 | 0.0125 | 0.0175 | 0.0056 | 0.0125 | 0.0071 | 0.0229 |
| 2 | house | 200 | 4 | 0.0331 | 20 | 0.0208 | 0.0046 | 0.0285 | 0.0378 | 0.0097 | 0.0250 | 0.0123 | 0.0539 |
| 2 | house | 200 | 5 | 0.9670 | 28 | 0.0262 | 0.0049 | 0.9620 | 0.9719 | 0.0101 | 0.9750 | 0.9408 | 0.9931 |
| 2 | house | 200 | 2 | 0.0000 | 0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
rt.RainPlot.E205.Block <- {
ggplot(data = avg.resp.long.E205.Block, aes(y = RespRate, x = as.factor(Resp), fill = Category)) +
geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .7) +
geom_point(aes(y = RespRate, color = Category), position = position_jitter(width = .15), size = .5, alpha = .8) +
geom_boxplot(aes(y = RespRate), width = .2, outlier.shape = NA, alpha = .7) +
geom_point(data = desc.resp.E205.Block, aes(y = Mean, x = Resp, color = Category), position = position_nudge(x = 0.3), size = 2.5) +
geom_errorbar(data = desc.resp.E205.Block, aes(y = Mean, ymin = Lower, ymax = Upper), position = position_nudge(x = 0.3), width = 0) +
facet_grid(Duration ~ blockID) +
# facet_grid(. ~ Type, switch = "x") + # create two panels to show data
# scale_colour_grey() + # start = .1, end = .6, color for the contour
# scale_fill_grey() + # start = .3, end = .6, color for the fill
# scale_color_brewer(palette = "Set1") + # palette = "RdBu"
# scale_fill_brewer(palette = "Set1") + # palette = "RdBu"
labs(title = "Reaction times for E205", x = "Stimulus Type", y = "Reaction times (ms)", fill = "Duration(ms)", color = "Duration(ms)") + # set the names for main, x and y axises and the legend
# coord_cartesian(ylim = c(0, 1.05)) + # set the limit for y axis
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
theme_bw() + # the backgroud color
plot_theme
}
rt.RainPlot.E205.Block
bloc.labels <- c("1st Half", "2nd Half")
names(bloc.labels) <- c(1, 2)
acc.ColuPlot.E205.Block <- {
ggplot(data = desc.resp.E205.Block, aes(y = Mean, x = Resp, fill = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Duration ~ blockID,
labeller = labeller(Duration = dura.labels,
blockID = bloc.labels)) +
geom_errorbar(mapping = aes(ymin = SE.lo, ymax = SE.hi), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
coord_cartesian(ylim = c(0,1.1)) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(title = "", x = "Responses", y = "Percentage", fill = "Category") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
erp_theme
}
acc.ColuPlot.E205.Block
df.ratio.E205.Block <- {
clean.beha.205 %>%
filter(Type == "intact") %>%
mutate(isFace = if_else(Category == "face", 1, 0)) %>%
group_by(SubjCode, blockID, Duration, Resp) %>%
summarize(ratio = mean(isFace), Count = n())
}
desc.ratio.E205.Block <- {
df.ratio.E205.Block %>%
group_by(blockID, Duration, Resp) %>%
summarize(Mean = mean(ratio)
, N = n()
, SD = if_else(N > 1, sd(ratio), 0)
, SE = SD/sqrt(N)
, SE.lo = Mean - SE, SE.hi = Mean + SE
, CI = SE * qt(0.975,N)
, Median = median(ratio)
, Lower = Mean - SD, Upper = Mean + SD # for rainclound plot
) %>%
ungroup()
}
ratio.ColuPlot.E205.Block <- {
ggplot(data = desc.ratio.E205.Block, aes(y = Mean, x = Resp)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_col(position = "dodge", color = "black", alpha = .7) + # position of columns and countour of columns
facet_grid(Duration ~ blockID,
labeller = labeller(Duration = dura.labels,
blockID = bloc.labels)) +
geom_errorbar(mapping = aes(ymin = SE.lo, ymax = SE.hi), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5,
position = position_dodge(width=0.9)) +
geom_hline(yintercept = c(0.5, 1), linetype = 5, alpha = 0.5) + # add the line for 0.5 and 1 (y)
coord_cartesian(ylim = c(0,1.1)) + # set the limit for y axis
scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(title = "", x = "Responses", y = "Percentage") + # set the names for main, x and y axises Subjective Responses for E205
scale_fill_grey() + # set the color for columns
# geom_text(label = c("***", "***", "", ""), color = c("red", "red", "red", "red"), size = 6, nudge_y = 0.05, nudge_x = rep(c(-0.225, 0.225), 2)) + # add starts to the significant columns
plot_theme
}
ratio.ColuPlot.E205.Block
df.erp.high.E2 <- {
df.erp.SC.E2 %>%
filter(Confidence == "high") %>%
mutate(DuraBlock = paste(Duration, Block_ST_table, sep = "_")) %>%
sdif_coding_E205_erp_conf()
}
# only keep the data for amplitudes of the P1 (correct and incorrect erps)
df.erp.P1.high.E2 <- {
df.erp.high.E2 %>%
filter(Component == "P1")
}
if (saveCSV) {
output.erp.P1.high.E2 <- file.path(folder.nesi, "E205_erp_P1_high.RData")
save(df.erp.P1.high.E2, file = output.erp.P1.high.E2)
}
# the structure of this dataset
head(df.erp.P1.high.E2, 10)
lmm.max.P1.high.E2 <- lmer(MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 + Hemi_C | Stimuli),
data = df.erp.P1.high.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e5)))
## start par. = 1 0 1 1 0 0 0 1 0 0 1 0 1 fn = 104824.6
## At return
## eval: 828 fn: 103991.60 par: 0.0607829 0.0132273 0.00000 0.368839 -0.0559292 0.0565168 0.0161722 0.107212 0.0588553 -0.00506718 0.178642 -0.0172574 1.00313e-07
## boundary (singular) fit: see ?isSingular
print(summary(lmm.max.P1.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category * Hemisphere * DuraBlock + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 + Hemi_C | Stimuli)
## Data: df.erp.P1.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 103991.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3462 -0.5802 -0.0117 0.5811 6.3206
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.059709 0.24436
## Hemi_C 0.002828 0.05318 1.00
## SubjCode (Intercept) 2.198646 1.48278
## Cate_C 0.236322 0.48613 -0.46
## Hemi_C 0.623365 0.78953 0.29 0.13
## Cate_Hemi 0.009455 0.09724 0.67 -0.49 -0.52
## Residual 16.161464 4.02013
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.379e+00 2.848e-01 2.794e+01 4.843 4.29e-05 ***
## Category2-1 -3.119e-03 1.354e-01 4.647e+01 -0.023 0.982
## Hemisphere2-1 1.101e-01 1.695e-01 3.103e+01 0.650 0.521
## DuraBlock2-1 -3.110e-02 9.951e-02 1.778e+04 -0.313 0.755
## DuraBlock3-2 1.062e-01 9.877e-02 1.193e+04 1.075 0.282
## Category2-1:Hemisphere2-1 -1.317e-01 1.491e-01 1.563e+03 -0.883 0.377
## Category2-1:DuraBlock2-1 -2.797e-02 1.988e-01 1.724e+04 -0.141 0.888
## Category2-1:DuraBlock3-2 1.071e-01 1.947e-01 5.630e+03 0.550 0.582
## Hemisphere2-1:DuraBlock2-1 4.922e-02 1.975e-01 1.801e+04 0.249 0.803
## Hemisphere2-1:DuraBlock3-2 1.643e-02 1.924e-01 1.163e+04 0.085 0.932
## Category2-1:Hemisphere2-1:DuraBlock2-1 -1.170e-02 3.940e-01 1.818e+04 -0.030 0.976
## Category2-1:Hemisphere2-1:DuraBlock3-2 -1.721e-01 3.737e-01 1.591e+04 -0.461 0.645
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
lmm.zcp.P1.high.E2 <- update(lmm.max.P1.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 + Hemi_C || Stimuli))
## start par. = 1 1 1 1 1 1 fn = 104824.6
## At return
## eval: 398 fn: 103999.59 par: 0.00000 0.0606616 3.15849e-07 0.196283 0.126126 0.371879
## boundary (singular) fit: see ?isSingular
print(summary(lmm.zcp.P1.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C +
## Cate_Hemi || SubjCode) + (1 + Hemi_C || Stimuli) + Category:Hemisphere +
## Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.P1.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 103999.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3295 -0.5808 -0.0128 0.5814 6.3179
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Hemi_C 0.00000 0.0000
## Stimuli.1 (Intercept) 0.05947 0.2439
## SubjCode Cate_Hemi 0.00000 0.0000
## SubjCode.1 Hemi_C 0.62267 0.7891
## SubjCode.2 Cate_C 0.25710 0.5071
## SubjCode.3 (Intercept) 2.23511 1.4950
## Residual 16.16202 4.0202
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.382e+00 2.872e-01 2.798e+01 4.811 4.66e-05 ***
## Category2-1 6.224e-03 1.392e-01 4.769e+01 0.045 0.965
## Hemisphere2-1 1.018e-01 1.696e-01 3.174e+01 0.601 0.552
## DuraBlock2-1 -3.233e-02 9.964e-02 1.821e+04 -0.324 0.746
## DuraBlock3-2 1.035e-01 9.925e-02 1.579e+04 1.043 0.297
## Category2-1:Hemisphere2-1 -1.448e-01 1.482e-01 1.785e+04 -0.977 0.328
## Category2-1:DuraBlock2-1 -2.736e-02 1.991e-01 1.794e+04 -0.137 0.891
## Category2-1:DuraBlock3-2 9.423e-02 1.962e-01 9.184e+03 0.480 0.631
## Hemisphere2-1:DuraBlock2-1 4.850e-02 1.976e-01 1.831e+04 0.245 0.806
## Hemisphere2-1:DuraBlock3-2 2.796e-02 1.929e-01 1.636e+04 0.145 0.885
## Category2-1:Hemisphere2-1:DuraBlock2-1 -6.056e-03 3.941e-01 1.832e+04 -0.015 0.988
## Category2-1:Hemisphere2-1:DuraBlock3-2 -1.563e-01 3.742e-01 1.830e+04 -0.418 0.676
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.zcp.P1.high.E2))
## $Stimuli
## Importance of components:
## [,1] [,2]
## Standard deviation 0.06066 0
## Proportion of Variance 1.00000 0
## Cumulative Proportion 1.00000 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.3719 0.1963 0.12613 0
## Proportion of Variance 0.7176 0.1999 0.08254 0
## Cumulative Proportion 0.7176 0.9175 1.00000 1
lmm.rdc.P1.high.E2 <- update(lmm.zcp.P1.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C || SubjCode) +
(1 | Stimuli))
## start par. = 1 1 1 1 fn = 104459.8
## At return
## eval: 149 fn: 103999.59 par: 0.0606616 0.196283 0.126126 0.371879
print(summary(lmm.rdc.P1.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C ||
## SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraBlock +
## Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.P1.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 103999.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3295 -0.5808 -0.0128 0.5814 6.3179
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli (Intercept) 0.05947 0.2439
## SubjCode Hemi_C 0.62267 0.7891
## SubjCode.1 Cate_C 0.25710 0.5070
## SubjCode.2 (Intercept) 2.23511 1.4950
## Residual 16.16202 4.0202
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.382e+00 2.872e-01 2.798e+01 4.811 4.66e-05 ***
## Category2-1 6.224e-03 1.392e-01 4.769e+01 0.045 0.965
## Hemisphere2-1 1.018e-01 1.696e-01 3.174e+01 0.601 0.552
## DuraBlock2-1 -3.233e-02 9.964e-02 1.821e+04 -0.324 0.746
## DuraBlock3-2 1.035e-01 9.925e-02 1.579e+04 1.043 0.297
## Category2-1:Hemisphere2-1 -1.448e-01 1.482e-01 1.785e+04 -0.977 0.328
## Category2-1:DuraBlock2-1 -2.736e-02 1.991e-01 1.794e+04 -0.137 0.891
## Category2-1:DuraBlock3-2 9.423e-02 1.962e-01 9.184e+03 0.480 0.631
## Hemisphere2-1:DuraBlock2-1 4.850e-02 1.976e-01 1.831e+04 0.245 0.806
## Hemisphere2-1:DuraBlock3-2 2.796e-02 1.929e-01 1.636e+04 0.145 0.885
## Category2-1:Hemisphere2-1:DuraBlock2-1 -6.056e-03 3.941e-01 1.832e+04 -0.015 0.988
## Category2-1:Hemisphere2-1:DuraBlock3-2 -1.563e-01 3.742e-01 1.830e+04 -0.418 0.676
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.P1.high.E2, lmm.zcp.P1.high.E2, refit = FALSE)
lmm.etd.P1.high.E2 <- update(lmm.rdc.P1.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C | SubjCode) +
(1 | Stimuli))
## start par. = 1 1 0 0 1 0 1 fn = 104459.8
## At return
## eval: 348 fn: 103992.90 par: 0.0606161 0.368821 -0.0560006 0.0540839 0.107121 0.0594284 0.180417
print(summary(lmm.etd.P1.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C |
## SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraBlock +
## Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.P1.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 103992.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3330 -0.5802 -0.0119 0.5800 6.3237
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.05939 0.2437
## SubjCode (Intercept) 2.19869 1.4828
## Cate_C 0.23616 0.4860 -0.46
## Hemi_C 0.63049 0.7940 0.27 0.14
## Residual 16.16335 4.0204
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.379e+00 2.848e-01 2.794e+01 4.842 4.29e-05 ***
## Category2-1 -3.143e-03 1.353e-01 4.641e+01 -0.023 0.982
## Hemisphere2-1 1.061e-01 1.703e-01 3.162e+01 0.623 0.538
## DuraBlock2-1 -3.146e-02 9.952e-02 1.778e+04 -0.316 0.752
## DuraBlock3-2 1.067e-01 9.878e-02 1.192e+04 1.081 0.280
## Category2-1:Hemisphere2-1 -1.444e-01 1.481e-01 1.772e+04 -0.975 0.330
## Category2-1:DuraBlock2-1 -2.871e-02 1.988e-01 1.724e+04 -0.144 0.885
## Category2-1:DuraBlock3-2 1.074e-01 1.947e-01 5.619e+03 0.552 0.581
## Hemisphere2-1:DuraBlock2-1 4.903e-02 1.976e-01 1.830e+04 0.248 0.804
## Hemisphere2-1:DuraBlock3-2 2.236e-02 1.928e-01 1.603e+04 0.116 0.908
## Category2-1:Hemisphere2-1:DuraBlock2-1 -9.847e-03 3.941e-01 1.832e+04 -0.025 0.980
## Category2-1:Hemisphere2-1:DuraBlock3-2 -1.557e-01 3.742e-01 1.829e+04 -0.416 0.677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.etd.P1.high.E2, lmm.rdc.P1.high.E2, refit = FALSE)
lmm.opt.P1.high.E205 <- lmm.etd.P1.high.E2
summary(lmm.opt.P1.high.E205)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C |
## SubjCode) + (1 | Stimuli) + Category:Hemisphere + Category:DuraBlock +
## Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.P1.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 103992.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3330 -0.5802 -0.0119 0.5800 6.3237
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.05939 0.2437
## SubjCode (Intercept) 2.19869 1.4828
## Cate_C 0.23616 0.4860 -0.46
## Hemi_C 0.63049 0.7940 0.27 0.14
## Residual 16.16335 4.0204
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.379e+00 2.848e-01 2.794e+01 4.842 4.29e-05 ***
## Category2-1 -3.143e-03 1.353e-01 4.641e+01 -0.023 0.982
## Hemisphere2-1 1.061e-01 1.703e-01 3.162e+01 0.623 0.538
## DuraBlock2-1 -3.146e-02 9.952e-02 1.778e+04 -0.316 0.752
## DuraBlock3-2 1.067e-01 9.878e-02 1.192e+04 1.081 0.280
## Category2-1:Hemisphere2-1 -1.444e-01 1.481e-01 1.772e+04 -0.975 0.330
## Category2-1:DuraBlock2-1 -2.871e-02 1.988e-01 1.724e+04 -0.144 0.885
## Category2-1:DuraBlock3-2 1.074e-01 1.947e-01 5.619e+03 0.552 0.581
## Hemisphere2-1:DuraBlock2-1 4.903e-02 1.976e-01 1.830e+04 0.248 0.804
## Hemisphere2-1:DuraBlock3-2 2.236e-02 1.928e-01 1.603e+04 0.116 0.908
## Category2-1:Hemisphere2-1:DuraBlock2-1 -9.847e-03 3.941e-01 1.832e+04 -0.025 0.980
## Category2-1:Hemisphere2-1:DuraBlock3-2 -1.557e-01 3.742e-01 1.829e+04 -0.416 0.677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ctg2-1 Hms2-1 DrB2-1 DrB3-2 Ct2-1:H2-1 C2-1:DB2 C2-1:DB3 H2-1:DB2 H2-1:DB3 C2-1:H2-1:DB2
## Category2-1 -0.267
## Hemisphr2-1 0.238 0.089
## DuraBlck2-1 0.036 0.027 -0.001
## DuraBlck3-2 -0.082 -0.193 -0.004 -0.660
## Ctg2-1:H2-1 0.001 0.005 0.184 -0.001 -0.004
## Ct2-1:DB2-1 0.007 0.158 -0.001 0.442 -0.252 -0.001
## Ct2-1:DB3-2 -0.047 -0.321 -0.005 -0.254 0.333 -0.004 -0.674
## Hm2-1:DB2-1 0.000 -0.001 0.131 0.001 0.001 0.070 0.001 0.001
## Hm2-1:DB3-2 -0.001 -0.007 -0.242 0.001 0.004 -0.279 0.001 0.007 -0.685
## C2-1:H2-1:DB2 0.000 0.000 0.033 0.001 0.000 0.316 0.000 0.000 0.439 -0.268
## C2-1:H2-1:DB3 -0.001 -0.002 -0.125 0.000 0.002 -0.468 0.000 0.002 -0.273 0.309 -0.714
anova(lmm.opt.P1.high.E205)
qqplot_lmer(lmm.opt.P1.high.E205)
emm.P1.high.E205 <- emmeans(lmm.opt.P1.high.E205, ~ DuraBlock | Category + Hemisphere)
emm.P1.high.E205
## Category = face, Hemisphere = Left:
## DuraBlock emmean SE df lower.CL upper.CL
## 17_1 1.320018 0.3099943 30.62 0.6874649 1.952572
## 17_2 1.275938 0.3261743 37.54 0.6153664 1.936510
## 200_2 1.278862 0.3116806 31.33 0.6434601 1.914265
##
## Category = house, Hemisphere = Left:
## DuraBlock emmean SE df lower.CL upper.CL
## 17_1 1.343170 0.2970638 47.27 0.7456435 1.940697
## 17_2 1.275302 0.3226546 66.22 0.6311419 1.919462
## 200_2 1.463480 0.2674223 31.68 0.9185459 2.008415
##
## Category = face, Hemisphere = Right:
## DuraBlock emmean SE df lower.CL upper.CL
## 17_1 1.428909 0.3422225 29.75 0.7297558 2.128062
## 17_2 1.438786 0.3569574 35.22 0.7142875 2.163284
## 200_2 1.541903 0.3438713 30.36 0.8399701 2.243835
##
## Category = house, Hemisphere = Right:
## DuraBlock emmean SE df lower.CL upper.CL
## 17_1 1.366141 0.3376954 41.64 0.6844687 2.047814
## 17_2 1.342383 0.3600515 53.99 0.6205191 2.064247
## 200_2 1.475083 0.3105991 30.11 0.8408562 2.109311
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
# emmip(emm.P1.high.E205, Category ~ DuraBlock | Hemisphere, CIs = TRUE)
durablock.labels = c("17ms\nhalf 1", "17ms\nhalf 2", "200ms\nhalf 2")
p1.high.LinePlot = {
ggplot(data = data.frame(emm.P1.high.E205), aes(y = emmean, x = DuraBlock, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit.p1) + # set the limit for y axis
scale_x_discrete(labels = durablock.labels) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Subjective Confidence", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# p1.high.LinePlot.slide <- {
# p1.high.LinePlot +
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('P1_high.png', p1.high.LinePlot.slide, width = 7.5, height = 4)
p1.high.LinePlot
# contrast(emmeans(lmm.opt.P1.high.E205, ~ DuraBlock | Category + Hemisphere), "pairwise", simple = "each", combine = TRUE)
contra.high.P1.gen <- contrast(emmeans(lmm.opt.P1.high.E205, ~ DuraBlock | Hemisphere + Category), #
"pairwise", adjust = "none")
contra.high.P1.gen[1:3]
## contrast Hemisphere Category estimate SE df t.ratio p.value
## 17_1 - 17_2 Left face 0.04407999 0.1482603 18312.03 0.297 0.7662
## 17_1 - 200_2 Left face 0.04115600 0.1227113 15818.16 0.335 0.7373
## 17_2 - 200_2 Left face -0.00292399 0.1591442 17711.71 -0.018 0.9853
##
## Degrees-of-freedom method: satterthwaite
confint(contra.high.P1.gen[1:3])
contra.high.P1.gen[4:6]
## contrast Hemisphere Category estimate SE df t.ratio p.value
## 17_1 - 17_2 Right face -0.0098768 0.1481876 18277.67 -0.067 0.9469
## 17_1 - 200_2 Right face -0.1129938 0.1223560 14516.06 -0.923 0.3558
## 17_2 - 200_2 Right face -0.1031170 0.1589035 17271.69 -0.649 0.5164
##
## Degrees-of-freedom method: satterthwaite
confint(contra.high.P1.gen[4:6])
contra.high.P1.gen.hemi <- contrast(emmeans(lmm.opt.P1.high.E205, ~ DuraBlock + Hemisphere | Category), #
interaction = "pairwise", adjust = "none")
contra.high.P1.gen.hemi
## Category = face:
## DuraBlock_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## 17_1 - 17_2 Left - Right 0.05395679 0.2091158 18274.58 0.258 0.7964
## 17_1 - 200_2 Left - Right 0.15414976 0.1710029 15420.21 0.901 0.3674
## 17_2 - 200_2 Left - Right 0.10019297 0.2233042 17648.80 0.449 0.6537
##
## Category = house:
## DuraBlock_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## 17_1 - 17_2 Left - Right 0.04410988 0.3347093 18314.99 0.132 0.8952
## 17_1 - 200_2 Left - Right -0.01136829 0.2472187 14349.79 -0.046 0.9633
## 17_2 - 200_2 Left - Right -0.05547817 0.3073727 17460.49 -0.180 0.8568
##
## Degrees-of-freedom method: satterthwaite
contra.high.P1.spec <- contrast(emmeans(lmm.opt.P1.high.E205, ~ Category + DuraBlock | Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni")
contra.high.P1.spec
## Hemisphere = Left:
## Category_pairwise DuraBlock_pairwise estimate SE df t.ratio p.value
## face - house 17_1 - 17_2 -0.02378829 0.2798066 18069.83 -0.085 1.0000
## face - house 17_1 - 200_2 0.16146633 0.2148371 3745.19 0.752 1.0000
## face - house 17_2 - 200_2 0.18525462 0.2697033 10769.68 0.687 1.0000
##
## Hemisphere = Right:
## Category_pairwise DuraBlock_pairwise estimate SE df t.ratio p.value
## face - house 17_1 - 17_2 -0.03363520 0.2799462 18184.19 -0.120 1.0000
## face - house 17_1 - 200_2 -0.00405172 0.2160355 4863.00 -0.019 1.0000
## face - house 17_2 - 200_2 0.02958348 0.2703208 12065.72 0.109 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 3 tests
confint(contra.high.P1.spec)
contra.high.P1.spec.hemi <- contrast(emmeans(lmm.opt.P1.high.E205, ~ Category + DuraBlock + Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni")
contra.high.P1.spec.hemi
## Category_pairwise DuraBlock_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## face - house 17_1 - 17_2 Left - Right 0.00984691 0.3940770 18315.17 0.025 1.0000
## face - house 17_1 - 200_2 Left - Right 0.16551805 0.2910815 18026.44 0.569 1.0000
## face - house 17_2 - 200_2 Left - Right 0.15567114 0.3741917 18291.34 0.416 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 3 tests
# only keep the data for amplitudes of the N170 (correct and incorrect erps)
df.erp.N170.high.E2 <- {
df.erp.high.E2 %>%
filter(Component == "N170")
}
if (saveCSV) {
output.erp.N170.high.E2 <- file.path(folder.nesi, "E205_erp_N170_high.RData")
save(df.erp.N170.high.E2, file = output.erp.N170.high.E2)
}
# the structure of this dataset
head(df.erp.N170.high.E2, 10)
lmm.max.N170.high.E2 <- lmer(MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 + Hemi_C | Stimuli),
data = df.erp.N170.high.E2,
verbose = TRUE,
control = lmerControl(optimizer = "bobyqa", # nloptwrap Nelder_Mead
optCtrl = list(maxfun = 1e5)))
## start par. = 1 0 1 1 0 0 0 1 0 0 1 0 1 fn = 105474.1
## At return
## eval: 1055 fn: 104763.30 par: 0.0573010 0.0277337 0.00000 0.454977 -0.134782 -0.0159484 -0.0974501 0.235464 -0.0364190 -0.107638 0.347869 -0.233303 0.243925
## boundary (singular) fit: see ?isSingular
print(summary(lmm.max.N170.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category * Hemisphere * DuraBlock + (1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) + (1 + Hemi_C | Stimuli)
## Data: df.erp.N170.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 104763.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6334 -0.5935 -0.0012 0.5937 6.1298
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.05502 0.2346
## Hemi_C 0.01289 0.1135 1.00
## SubjCode (Intercept) 3.46891 1.8625
## Cate_C 1.23353 1.1106 -0.50
## Hemi_C 2.05438 1.4333 -0.05 -0.07
## Cate_Hemi 2.26249 1.5042 -0.27 -0.12 -0.59
## Residual 16.75767 4.0936
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.319e+00 3.558e-01 2.775e+01 -6.517 4.80e-07 ***
## Category2-1 1.625e+00 2.339e-01 3.373e+01 6.949 5.42e-08 ***
## Hemisphere2-1 -3.835e-01 2.846e-01 2.969e+01 -1.348 0.188
## DuraBlock2-1 -8.716e-01 1.016e-01 1.831e+04 -8.581 < 2e-16 ***
## DuraBlock3-2 1.019e+00 1.014e-01 1.710e+04 10.052 < 2e-16 ***
## Category2-1:Hemisphere2-1 4.113e-01 3.317e-01 3.519e+01 1.240 0.223
## Category2-1:DuraBlock2-1 -1.442e-01 2.030e-01 1.827e+04 -0.710 0.478
## Category2-1:DuraBlock3-2 1.909e+00 2.018e-01 1.457e+04 9.460 < 2e-16 ***
## Hemisphere2-1:DuraBlock2-1 -1.955e-01 2.024e-01 1.765e+04 -0.966 0.334
## Hemisphere2-1:DuraBlock3-2 4.096e-02 1.997e-01 1.025e+04 0.205 0.837
## Category2-1:Hemisphere2-1:DuraBlock2-1 4.030e-01 4.045e-01 1.739e+04 0.996 0.319
## Category2-1:Hemisphere2-1:DuraBlock3-2 3.517e-03 3.962e-01 6.953e+03 0.009 0.993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
lmm.zcp.N170.high.E2 <- update(lmm.max.N170.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 + Hemi_C || Stimuli))
## start par. = 1 1 1 1 1 1 fn = 105474.1
## At return
## eval: 483 fn: 104787.06 par: 1.26220e-07 0.0562734 0.374364 0.356540 0.272915 0.456967
## boundary (singular) fit: see ?isSingular
print(summary(lmm.zcp.N170.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C +
## Cate_Hemi || SubjCode) + (1 + Hemi_C || Stimuli) + Category:Hemisphere +
## Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.N170.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 104787.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6562 -0.5916 0.0000 0.5928 6.1351
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli Hemi_C 0.00000 0.0000
## Stimuli.1 (Intercept) 0.05308 0.2304
## SubjCode Cate_Hemi 2.34918 1.5327
## SubjCode.1 Hemi_C 2.13081 1.4597
## SubjCode.2 Cate_C 1.24849 1.1174
## SubjCode.3 (Intercept) 3.50023 1.8709
## Residual 16.76210 4.0942
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.316e+00 3.574e-01 2.775e+01 -6.480 5.30e-07 ***
## Category2-1 1.639e+00 2.352e-01 3.368e+01 6.968 5.16e-08 ***
## Hemisphere2-1 -3.760e-01 2.898e-01 2.981e+01 -1.297 0.205
## DuraBlock2-1 -8.721e-01 1.017e-01 1.836e+04 -8.579 < 2e-16 ***
## DuraBlock3-2 1.015e+00 1.017e-01 1.820e+04 9.981 < 2e-16 ***
## Category2-1:Hemisphere2-1 4.367e-01 3.389e-01 3.544e+01 1.289 0.206
## Category2-1:DuraBlock2-1 -1.426e-01 2.032e-01 1.835e+04 -0.702 0.483
## Category2-1:DuraBlock3-2 1.888e+00 2.028e-01 1.705e+04 9.314 < 2e-16 ***
## Hemisphere2-1:DuraBlock2-1 -1.984e-01 2.029e-01 1.829e+04 -0.978 0.328
## Hemisphere2-1:DuraBlock3-2 3.172e-02 2.019e-01 1.734e+04 0.157 0.875
## Category2-1:Hemisphere2-1:DuraBlock2-1 4.145e-01 4.057e-01 1.823e+04 1.022 0.307
## Category2-1:Hemisphere2-1:DuraBlock3-2 -4.034e-02 4.019e-01 1.438e+04 -0.100 0.920
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
summary(rePCA(lmm.zcp.N170.high.E2))
## $Stimuli
## Importance of components:
## [,1] [,2]
## Standard deviation 0.05627 0
## Proportion of Variance 1.00000 0
## Cumulative Proportion 1.00000 1
##
## $SubjCode
## Importance of components:
## [,1] [,2] [,3] [,4]
## Standard deviation 0.4570 0.3744 0.3565 0.2729
## Proportion of Variance 0.3793 0.2545 0.2309 0.1353
## Cumulative Proportion 0.3793 0.6338 0.8647 1.0000
lmm.rdc.N170.high.E2 <- update(lmm.zcp.N170.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi || SubjCode) +
(1 | Stimuli))
## start par. = 1 1 1 1 1 fn = 105207.8
## At return
## eval: 259 fn: 104787.06 par: 0.0562734 0.374363 0.356541 0.272917 0.456967
print(summary(lmm.rdc.N170.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C +
## Cate_Hemi || SubjCode) + (1 | Stimuli) + Category:Hemisphere +
## Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.N170.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 104787.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6562 -0.5916 0.0000 0.5928 6.1351
##
## Random effects:
## Groups Name Variance Std.Dev.
## Stimuli (Intercept) 0.05308 0.2304
## SubjCode Cate_Hemi 2.34917 1.5327
## SubjCode.1 Hemi_C 2.13082 1.4597
## SubjCode.2 Cate_C 1.24850 1.1174
## SubjCode.3 (Intercept) 3.50024 1.8709
## Residual 16.76210 4.0942
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.316e+00 3.574e-01 2.775e+01 -6.480 5.30e-07 ***
## Category2-1 1.639e+00 2.352e-01 3.368e+01 6.968 5.16e-08 ***
## Hemisphere2-1 -3.760e-01 2.898e-01 2.981e+01 -1.297 0.205
## DuraBlock2-1 -8.721e-01 1.017e-01 1.836e+04 -8.579 < 2e-16 ***
## DuraBlock3-2 1.015e+00 1.017e-01 1.820e+04 9.981 < 2e-16 ***
## Category2-1:Hemisphere2-1 4.367e-01 3.389e-01 3.544e+01 1.289 0.206
## Category2-1:DuraBlock2-1 -1.426e-01 2.032e-01 1.835e+04 -0.702 0.483
## Category2-1:DuraBlock3-2 1.888e+00 2.028e-01 1.705e+04 9.314 < 2e-16 ***
## Hemisphere2-1:DuraBlock2-1 -1.984e-01 2.029e-01 1.829e+04 -0.978 0.328
## Hemisphere2-1:DuraBlock3-2 3.172e-02 2.019e-01 1.734e+04 0.157 0.875
## Category2-1:Hemisphere2-1:DuraBlock2-1 4.145e-01 4.057e-01 1.823e+04 1.022 0.307
## Category2-1:Hemisphere2-1:DuraBlock3-2 -4.034e-02 4.019e-01 1.438e+04 -0.100 0.920
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.rdc.N170.high.E2, lmm.zcp.N170.high.E2, refit = FALSE)
lmm.etd.N170.high.E2 <- update(lmm.rdc.N170.high.E2,
formula = MeanAmp ~ Category * Hemisphere * DuraBlock +
(1 + Cate_C + Hemi_C + Cate_Hemi | SubjCode) +
(1 | Stimuli))
## start par. = 1 1 0 0 0 1 0 0 1 0 1 fn = 105207.8
## At return
## eval: 903 fn: 104765.64 par: 0.0562480 0.454957 -0.134602 -0.0160593 -0.0976551 0.235449 -0.0361790 -0.107527 0.347685 -0.233521 0.243953
print(summary(lmm.etd.N170.high.E2), corr = FALSE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C +
## Cate_Hemi | SubjCode) + (1 | Stimuli) + Category:Hemisphere +
## Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.N170.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 104765.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6367 -0.5923 -0.0010 0.5935 6.1423
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.05303 0.2303
## SubjCode (Intercept) 3.46964 1.8627
## Cate_C 1.23296 1.1104 -0.50
## Hemi_C 2.05262 1.4327 -0.05 -0.07
## Cate_Hemi 2.26537 1.5051 -0.27 -0.12 -0.59
## Residual 16.76265 4.0942
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.319e+00 3.558e-01 2.774e+01 -6.517 4.82e-07 ***
## Category2-1 1.626e+00 2.336e-01 3.361e+01 6.960 5.34e-08 ***
## Hemisphere2-1 -3.885e-01 2.842e-01 2.957e+01 -1.367 0.182
## DuraBlock2-1 -8.726e-01 1.016e-01 1.832e+04 -8.590 < 2e-16 ***
## DuraBlock3-2 1.019e+00 1.014e-01 1.712e+04 10.054 < 2e-16 ***
## Category2-1:Hemisphere2-1 4.041e-01 3.309e-01 3.480e+01 1.221 0.230
## Category2-1:DuraBlock2-1 -1.453e-01 2.031e-01 1.828e+04 -0.716 0.474
## Category2-1:DuraBlock3-2 1.908e+00 2.018e-01 1.459e+04 9.454 < 2e-16 ***
## Hemisphere2-1:DuraBlock2-1 -1.971e-01 2.024e-01 1.765e+04 -0.974 0.330
## Hemisphere2-1:DuraBlock3-2 5.056e-02 1.997e-01 1.025e+04 0.253 0.800
## Category2-1:Hemisphere2-1:DuraBlock2-1 3.992e-01 4.045e-01 1.739e+04 0.987 0.324
## Category2-1:Hemisphere2-1:DuraBlock3-2 1.533e-02 3.962e-01 6.960e+03 0.039 0.969
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm.etd.N170.high.E2, lmm.rdc.N170.high.E2, refit = FALSE)
lmm.opt.N170.high.E205 <- lmm.etd.N170.high.E2
summary(lmm.opt.N170.high.E205)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: MeanAmp ~ Category + Hemisphere + DuraBlock + (1 + Cate_C + Hemi_C +
## Cate_Hemi | SubjCode) + (1 | Stimuli) + Category:Hemisphere +
## Category:DuraBlock + Hemisphere:DuraBlock + Category:Hemisphere:DuraBlock
## Data: df.erp.N170.high.E2
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
##
## REML criterion at convergence: 104765.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.6367 -0.5923 -0.0010 0.5935 6.1423
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Stimuli (Intercept) 0.05303 0.2303
## SubjCode (Intercept) 3.46964 1.8627
## Cate_C 1.23296 1.1104 -0.50
## Hemi_C 2.05262 1.4327 -0.05 -0.07
## Cate_Hemi 2.26537 1.5051 -0.27 -0.12 -0.59
## Residual 16.76265 4.0942
## Number of obs: 18456, groups: Stimuli, 80; SubjCode, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.319e+00 3.558e-01 2.774e+01 -6.517 4.82e-07 ***
## Category2-1 1.626e+00 2.336e-01 3.361e+01 6.960 5.34e-08 ***
## Hemisphere2-1 -3.885e-01 2.842e-01 2.957e+01 -1.367 0.182
## DuraBlock2-1 -8.726e-01 1.016e-01 1.832e+04 -8.590 < 2e-16 ***
## DuraBlock3-2 1.019e+00 1.014e-01 1.712e+04 10.054 < 2e-16 ***
## Category2-1:Hemisphere2-1 4.041e-01 3.309e-01 3.480e+01 1.221 0.230
## Category2-1:DuraBlock2-1 -1.453e-01 2.031e-01 1.828e+04 -0.716 0.474
## Category2-1:DuraBlock3-2 1.908e+00 2.018e-01 1.459e+04 9.454 < 2e-16 ***
## Hemisphere2-1:DuraBlock2-1 -1.971e-01 2.024e-01 1.765e+04 -0.974 0.330
## Hemisphere2-1:DuraBlock3-2 5.056e-02 1.997e-01 1.025e+04 0.253 0.800
## Category2-1:Hemisphere2-1:DuraBlock2-1 3.992e-01 4.045e-01 1.739e+04 0.987 0.324
## Category2-1:Hemisphere2-1:DuraBlock3-2 1.533e-02 3.962e-01 6.960e+03 0.039 0.969
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ctg2-1 Hms2-1 DrB2-1 DrB3-2 Ct2-1:H2-1 C2-1:DB2 C2-1:DB3 H2-1:DB2 H2-1:DB3 C2-1:H2-1:DB2
## Category2-1 -0.418
## Hemisphr2-1 -0.044 -0.061
## DuraBlck2-1 0.029 0.013 0.001
## DuraBlck3-2 -0.070 -0.126 0.004 -0.651
## Ctg2-1:H2-1 -0.228 -0.101 -0.415 0.002 0.009
## Ct2-1:DB2-1 0.004 0.089 0.001 0.444 -0.245 0.002
## Ct2-1:DB3-2 -0.041 -0.207 0.005 -0.246 0.345 0.011 -0.656
## Hm2-1:DB2-1 0.000 0.001 0.075 -0.002 -0.001 0.024 -0.002 -0.001
## Hm2-1:DB3-2 0.002 0.006 -0.162 -0.001 -0.009 -0.154 -0.002 -0.010 -0.667
## C2-1:H2-1:DB2 0.000 0.002 0.014 -0.002 -0.002 0.131 -0.003 -0.002 0.440 -0.254
## C2-1:H2-1:DB3 0.002 0.008 -0.090 -0.001 -0.010 -0.267 -0.002 -0.013 -0.255 0.325 -0.674
anova(lmm.opt.N170.high.E205)
qqplot_lmer(lmm.opt.N170.high.E205)
emm.N170.high.E205 <- emmeans(lmm.opt.N170.high.E205, ~ DuraBlock | Category + Hemisphere)
emm.N170.high.E205
## Category = face, Hemisphere = Left:
## DuraBlock emmean SE df lower.CL upper.CL
## 17_1 -2.449953 0.4543250 28.82 -3.379405 -1.5205000
## 17_2 -3.051512 0.4658692 31.86 -4.000617 -2.1024077
## 200_2 -3.007625 0.4549179 28.98 -3.938059 -2.0771908
##
## Category = house, Hemisphere = Left:
## DuraBlock emmean SE df lower.CL upper.CL
## 17_1 -1.429557 0.4056330 38.52 -2.250358 -0.6087573
## 17_2 -2.376014 0.4231693 45.76 -3.227932 -1.5240962
## 200_2 -0.431652 0.3761813 28.76 -1.201308 0.3380050
##
## Category = face, Hemisphere = Right:
## DuraBlock emmean SE df lower.CL upper.CL
## 17_1 -2.790329 0.4736852 28.54 -3.759797 -1.8208619
## 17_2 -3.788595 0.4847723 31.31 -4.776894 -2.8002960
## 200_2 -3.701812 0.4743015 28.70 -4.672301 -2.7313228
##
## Category = house, Hemisphere = Right:
## DuraBlock emmean SE df lower.CL upper.CL
## 17_1 -1.637055 0.3452892 42.41 -2.333678 -0.9404319
## 17_2 -2.581063 0.3668471 54.40 -3.316423 -1.8457034
## 200_2 -0.578478 0.3141882 29.66 -1.220446 0.0634906
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
# emmip(emm.N170.high.E205, Category ~ DuraBlock | Hemisphere, CIs = TRUE) +
# scale_y_continuous(trans = "reverse")
n170.high.LinePlot = {
ggplot(data = data.frame(emm.N170.high.E205), aes(y = emmean, x = DuraBlock, color = Category, group = Category)) + # set the data, varialbes for x and y axises, and the variable for dividing data
geom_point(size = 2) +
geom_line(aes(linetype = Category), size = .85) + # position = "dodge", alpha = .7
scale_linetype_manual(values=c("solid", "dashed")) +
geom_errorbar(mapping = aes(ymin = lower.CL, ymax = upper.CL), linetype = 1, # set the error bar
show.legend = FALSE, width = 0.25, alpha = .5) + # , position = position_dodge(width=0.9)
facet_grid(. ~ Hemisphere, # switch = "x",
labeller = labeller(Hemisphere = hemiLabel)) +
scale_y_reverse(limits = ylimit.n170, breaks = seq(1, -5, -1)) + # set the limit for y axis
scale_x_discrete(labels = durablock.labels) +
# scale_y_continuous(expand= c(0, 0)) + # remove the space between columns and x axis
labs(x = "Duration and Subjective Confidence", y = expression(paste("Amplitudes (", mu, "V)")), fill = "Category") + # set the names for main, x and y axises title = "",
# scale_x_discrete(labels = c("left", "right")) +
# geom_text(label = c("", "", "", "", "***", "", "", ""), color = "red", size = 6, nudge_y = 0.5, nudge_x = 0.5) + # add starts to the significant columns
theme_bw() +
erp_theme
}
# n170.high.LinePlot.slide <- {
# n170.high.LinePlot +
# scale_y_reverse(limits = c(0.8, -4.8), breaks = seq(0, -4, -1)) + # set the limit for y axis
# theme(legend.position="right",
# axis.line = element_line(color="black", size = 3))
# }
# ggsave('N170_high.png', n170.high.LinePlot.slide, width = 7.5, height = 4)
n170.high.LinePlot
# contra.high.N170 <- contrast(emmeans(lmm.opt.N170.high.E205, ~ Category + DuraBlock), "pairwise", simple = "each")
#
# contra.high.N170
#
# confint(contra.high.N170)
contra.high.N170.gen <- contrast(emmeans(lmm.opt.N170.high.E205, ~ DuraBlock | Hemisphere + Category), #
"pairwise", adjust = "Bonferroni")
contra.high.N170.gen[1:3]
## contrast Hemisphere Category estimate SE df t.ratio p.value
## 17_1 - 17_2 Left face 0.6015598 0.1514381 18306.67 3.972 0.0002
## 17_1 - 200_2 Left face 0.5576724 0.1272494 17686.67 4.383 <.0001
## 17_2 - 200_2 Left face -0.0438874 0.1636515 18181.18 -0.268 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 3 tests
confint(contra.high.N170.gen[1:3])
contra.high.N170.gen[4:6]
## contrast Hemisphere Category estimate SE df t.ratio p.value
## 17_1 - 17_2 Right face 0.9982659 0.1514036 18303.96 6.593 <.0001
## 17_1 - 200_2 Right face 0.9114827 0.1270620 17334.05 7.174 <.0001
## 17_2 - 200_2 Right face -0.0867832 0.1635189 18079.17 -0.531 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 3 tests
confint(contra.high.N170.gen[4:6])
contra.high.N170.gen.hemi <- contrast(emmeans(lmm.opt.N170.high.E205, ~ DuraBlock + Hemisphere | Category), #
interaction = "pairwise", adjust = "Bonferroni")
contra.high.N170.gen.hemi
## Category = face:
## DuraBlock_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## 17_1 - 17_2 Left - Right -0.3967060 0.2140373 18283.98 -1.853 0.1915
## 17_1 - 200_2 Left - Right -0.3538102 0.1794494 16921.73 -1.972 0.1460
## 17_2 - 200_2 Left - Right 0.0428958 0.2310594 17949.51 0.186 1.0000
##
## Category = house:
## DuraBlock_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## 17_1 - 17_2 Left - Right 0.0024481 0.3434039 16785.99 0.007 1.0000
## 17_1 - 200_2 Left - Right 0.0606714 0.2719402 1429.15 0.223 1.0000
## 17_2 - 200_2 Left - Right 0.0582233 0.3237972 5134.82 0.180 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 3 tests
contra.high.N170.spec <- contrast(emmeans(lmm.opt.N170.high.E205, ~ Category + DuraBlock | Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni")
contra.high.N170.spec[1:3]
## Category_pairwise DuraBlock_pairwise Hemisphere estimate SE df t.ratio p.value
## face - house 17_1 - 17_2 Left -0.3448969 0.2869752 18209.03 -1.202 0.6883
## face - house 17_1 - 200_2 Left 1.5555782 0.2363195 7112.00 6.583 <.0001
## face - house 17_2 - 200_2 Left 1.9004751 0.2846656 13731.61 6.676 <.0001
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 3 tests
confint(contra.high.N170.spec[1:3])
contra.high.N170.spec[4:6]
## Category_pairwise DuraBlock_pairwise Hemisphere estimate SE df t.ratio p.value
## face - house 17_1 - 17_2 Right 0.0542572 0.2862422 17571.01 0.190 1.0000
## face - house 17_1 - 200_2 Right 1.9700598 0.2297485 2326.23 8.575 <.0001
## face - house 17_2 - 200_2 Right 1.9158026 0.2809292 7615.60 6.820 <.0001
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 3 tests
confint(contra.high.N170.spec[4:6])
contra.high.N170.spec.hemi <- contrast(emmeans(lmm.opt.N170.high.E205, ~ Category + DuraBlock + Hemisphere), #
interaction = "pairwise", adjust = "Bonferroni")
contra.high.N170.spec.hemi
## Category_pairwise DuraBlock_pairwise Hemisphere_pairwise estimate SE df t.ratio p.value
## face - house 17_1 - 17_2 Left - Right -0.3991541 0.4044924 17388.02 -0.987 0.9713
## face - house 17_1 - 200_2 Left - Right -0.4144816 0.3231755 2014.65 -1.283 0.5994
## face - house 17_2 - 200_2 Left - Right -0.0153275 0.3961737 6959.66 -0.039 1.0000
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: bonferroni method for 3 tests
For the analysis for every dependent variable, an optimal model was obtained by the following general steps:
step(fit, fixed.reduce = FALSE) from library(lmerTest) is used for this step. The algorithm could be found here.More details about this whole process chould be found here: Bates, D., Kliegl, R., Vasishth, S., & Baayen, H. (2015). Parsimonious Mixed Models. Retrieved from http://arxiv.org/abs/1506.04967.
stepThe outline of the algorithm of step(fit, reduce.fixed = false) function is:
Simplification of the random effects structure:
Let M be the linear mixed effects model specified by a user.
If there are random effects in M then go to 3, otherwise stop.
- For each random effect ri in M do:
Create a reduced model Mi by eliminating ri from M.
Calculate pi, the p value from the likelihood ratio test of comparing M to Mi.
Save pi and Mi.
Find pmax; the maximum of all pi and let Mmax denote the corresponding model.
Set M to Mmax. If pmax is guesser than α level then go back to 3, otherwise stop.
More deatils could be found: Kuznetsova, A., Brockhoff, P. B., & Christensen, R. H. B. (2017). lmerTest Package: Tests in Linear Mixed Effects Models. Journal of Statistical Software, 82(13), 1-26. http://doi.org/10.18637/jss.v082.i13 Page 8
# rstudioapi::versionInfo()
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggpubr_0.2.5 emmeans_1.4.2 lmerTest_3.1-0 lme4_1.1-21 Matrix_1.2-17 magrittr_1.5 forcats_0.4.0
## [8] stringr_1.4.0 dplyr_0.8.3 purrr_0.3.3 readr_1.3.1 tidyr_1.0.0 tibble_2.1.3 ggplot2_3.2.1
## [15] tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.3 lubridate_1.7.4 mvtnorm_1.0-11 lattice_0.20-38 assertthat_0.2.1 zeallot_0.1.0
## [7] digest_0.6.22 R6_2.4.0 cellranger_1.1.0 plyr_1.8.4 backports_1.1.5 evaluate_0.14
## [13] coda_0.19-3 highr_0.8 httr_1.4.1 pillar_1.4.2 rlang_0.4.1 lazyeval_0.2.2
## [19] readxl_1.3.1 rstudioapi_0.10 minqa_1.2.4 nloptr_1.2.1 rmarkdown_1.16 labeling_0.3
## [25] splines_3.6.1 munsell_0.5.0 broom_0.5.2 compiler_3.6.1 numDeriv_2016.8-1.1 modelr_0.1.5
## [31] xfun_0.10 pkgconfig_2.0.3 htmltools_0.4.0 tidyselect_0.2.5 crayon_1.3.4 withr_2.1.2
## [37] MASS_7.3-51.4 grid_3.6.1 nlme_3.1-140 jsonlite_1.6 xtable_1.8-4 gtable_0.3.0
## [43] lifecycle_0.1.0 scales_1.0.0 estimability_1.3 cli_1.1.0 stringi_1.4.3 ggsignif_0.6.0
## [49] reshape2_1.4.3 xml2_1.2.2 generics_0.0.2 vctrs_0.2.0 cowplot_1.0.0 boot_1.3-22
## [55] tools_3.6.1 glue_1.3.1 hms_0.5.2 yaml_2.2.0 colorspace_1.4-1 rvest_0.3.5
## [61] knitr_1.25 haven_2.2.0
A work by Haiyang Jin